Changeset 2888


Ignore:
Timestamp:
Feb 1, 2023, 8:19:12 PM (22 months ago)
Author:
llange
Message:

PEM

  • Following r-2886-5, debugging some issues ( conservation of H2O, water_reservoir not initialized, info_pem which was not working)
  • Cleaning of some redundant routines (e.g., stopping criterion of the PEM) and variables renamed (e.g., alpha_criterion now transofrmed in ice_criterion)
  • Ice table subroutine is now in a module and recalled "compute_ice_table_equilibrium" (v.s. dynamic ice table, efforts ongoing)
  • Abort_pem introduced to avoid just stopping the pem with raw "stop" in the codes
  • conf_pem has been improved so that values are not hard coded (e.g., geothermal flux, mass of water_reservoir)
  • Regolith thermal properties updated in a cleaner way

(Not atomic commit, sorry)
LL & RV

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 added
3 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90

    r2863 r2888  
    11  module adsorption_mod 
    22  implicit none
     3  LOGICAL adsorption_pem ! True by default, to compute adsorption/desorption. Read in  pem.def
     4  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
     5  real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
    36  contains
    47
    … …  
    1013!!!
    1114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     15
     16
     17  subroutine ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
     18
     19  implicit none
     20  integer,intent(in) :: ngrid ! number of atmospheric columns
     21  integer,intent(in) :: nslope ! number of slope within a mesh
     22  integer,intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
     23    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
     24    allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
     25  end subroutine ini_adsorption_h_PEM
     26
     27  subroutine end_adsorption_h_PEM
     28
     29  implicit none
     30    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
     31    if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
     32  end subroutine end_adsorption_h_PEM
     33
     34
     35
     36
     37
    1238
    1339  subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &

  • trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90

    r2855 r2888  
    3131
    3232         enddo
    33 
     33         if (ig.eq.400) then
     34         write(*,*) 'ig,islope,diff',ig,islope,diff_rho(:)
     35         write(*,*) 'rho surf',rhowatersurf_ave(ig,islope)
     36         write(*,*) 'rho soil',rhowatersoil_ave(ig,:,islope)
     37         endif
    3438         if(diff_rho(1) > 0) then                                   ! ice is at the surface
    3539           ice_table(ig,islope) = 0.

  • trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90

    r2885 r2888  
    1515  real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients [SI]
    1616  real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients [SI]
    17   real,save,allocatable :: alph_PEM(:,:,:)      ! (FC) alpha_k coefficients [SI]
    18   real,save,allocatable :: beta_PEM(:,:,:)      ! beta_k coefficients [SI]
     17  real,save,allocatable :: alph_PEM(:,:)      ! (FC) alpha_k coefficients [SI]
     18  real,save,allocatable :: beta_PEM(:,:)      ! beta_k coefficients [SI]
    1919  real,save :: mu_PEM                           ! mu coefficient [SI]
    20   real,parameter :: fluxgeo = 30e-3             ! Geothermal flux [W/m^2]
    21   real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
    22   real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
     20  real,save :: fluxgeo                               ! Geothermal flux [W/m^2]
    2321  LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
    2422  real,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
    25 
     23  real,save :: water_reservoir_nom
    2624contains
    2725
    … …  
    4038    allocate(coefq_PEM(0:nsoilmx_PEM-1))
    4139    allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
    42     allocate(alph_PEM(ngrid,nsoilmx_PEM-1,nslope))
    43     allocate(beta_PEM(ngrid,nsoilmx_PEM-1,nslope))
     40    allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
     41    allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
    4442    allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
    45     allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
    46     allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
     43    allocate(water_reservoir(ngrid))
    4744    allocate(water_reservoir(ngrid))
    4845  end subroutine ini_comsoil_h_PEM
    … …  
    6461    if (allocated(beta_PEM)) deallocate(beta_PEM)
    6562    if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
    66     if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
    67     if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
     63    if (allocated(water_reservoir)) deallocate(water_reservoir)
    6864    if (allocated(water_reservoir)) deallocate(water_reservoir)
    6965  end subroutine end_comsoil_h_PEM

  • trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90

    r2864 r2888  
    1414#endif
    1515 
    16   USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, alpha_criterion, &
     16  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, ice_criterion, ps_criterion, &
    1717                Max_iter_pem, evol_orbit_pem
    18   USE comsoil_h_pem, only: soil_pem
    19 
    20  
     18  USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom
     19  USE adsorption_mod,only: adsorption_pem
     20  CHARACTER(len=20),parameter :: modname ='conf_pem'
    2121!PEM parameter
    2222  year_bp_ini=0.
    … …  
    2626  CALL getin('dt_pem', dt_pem)
    2727
    28   alpha_criterion=0.2
    29   CALL getin('alpha_criterion', alpha_criterion)
     28  ice_criterion=0.2
     29  CALL getin('ice_criterion', ice_criterion)
     30
     31  ps_criterion = 0.15
     32  CALL getin('ps_criterion',ps_criterion)
    3033
    3134  evol_orbit_pem=.false.
    … …  
    3841  CALL getin('soil_pem', soil_pem)
    3942
     43  adsorption_pem = .true.
     44  CALL getin('adsorption_pem',adsorption_pem)
     45
     46  fluxgeo = 0.
     47  CALL getin('Fluxgeo_PEM',fluxgeo)
     48  print*,'Flux Geothermal is set to',fluxgeo
     49
     50  if ((not(soil_pem)).and.adsorption_pem) then
     51       print*,'Adsorption must be used when soil_pem = T'
     52       call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
     53        endif
     54 
     55
     56  if ((not(soil_pem)).and.fluxgeo.gt.0.) then
     57       print*,'Soil is not activated but Flux Geo > 0.'
     58       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.'",1)
     59        endif
     60 
     61
     62   water_reservoir_nom = 1e4
     63  CALL getin('water_reservoir_nom',water_reservoir_nom)
    4064  END SUBROUTINE conf_pem
    4165

  • trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90

    r2863 r2888  
    77  USE temps_mod_evol, ONLY: dt_pem
    88  use comslope_mod, ONLY: subslope_dist
    9 
     9  use criterion_pem_stop_mod, only: criterion_waterice_stop
    1010      IMPLICIT NONE
    1111
    … …  
    8080    print *, "Tendencies on ice increasing=", pos_tend
    8181    print *, "This can be due to the absence of water ice in the PCM run!!"
    82       call criterion_ice_stop_water_slope(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
     82      call criterion_waterice_stop(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
    8383      do i=1,ngrid
    8484         do islope=1,nslope

  • trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90

    r2886 r2888  
    11SUBROUTINE info_run_PEM(year_iter, criterion_stop)
    22
    3       IMPLICIT NONE
     3   IMPLICIT NONE
    44
    55  INTEGER, intent(in) :: year_iter, criterion_stop                  ! # of year and reason to stop
    66
    7  open(1,file='info_run_PEM.txt', status='new')
    8 
    9  write(1,*)  year_iter, criterion_stop
    10 
    11  close(1)
     7  logical :: exist
    128
    139
     10  inquire(file='info_run_PEM.txt', exist=exist)
     11  if (exist) then
     12    open(12, file='info_run_PEM.txt', status="old", position="append", action="write")
     13  else
     14    open(12, file='info_run_PEM.txt', status="new", action="write")
     15  end if
     16  write(12, *) year_iter, criterion_stop
     17  close(12)
     18
    1419END SUBROUTINE info_run_PEM

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

    r2885 r2888  
    8484      USE geometry_mod, only: latitude_deg
    8585      use conf_pem_mod, only: conf_pem
    86       use pemredem, only:  pemdem1
     86      use pemredem, only:  pemdem0,pemdem1
    8787      use co2glaciers_mod,only: co2glaciers_evol
     88      use criterion_pem_stop_mod,only: criterion_waterice_stop,criterion_co2_stop
    8889!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
    8990      use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
    90                               TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia         
     91                              TI_PEM,inertiedat_PEM,         & ! soil thermal inertia         
    9192                              tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
    92                               fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys,  & ! geothermal flux, mass of co2 & h2o in the regolith
    93                               water_reservoir
    94       use adsorption_mod, only : regolith_adsorption
     93                              fluxgeo,water_reservoir                             ! geothermal flux
     94      use adsorption_mod, only : regolith_adsorption,adsorption_pem, &   ! bool to check if adsorption, main subroutine
     95                                 ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays                               
     96                                 co2_adsorbded_phys, h2o_adsorbded_phys  ! mass of co2 and h2O adsorbded
    9597
    9698!!! For orbit parameters
    … …  
    98100      use orbit_param_criterion_mod, only : orbit_param_criterion
    99101      use recomp_orb_param_mod, only: recomp_orb_param
     102      use ice_table_mod, only: computeice_table_equilibrium
    100103
    101104
    … …  
    189192      LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
    190193      LOGICAL :: STOPPING_pressure
    191       INTEGER :: criterion_stop            !Which criterion is reached : 1=surf h20, 2=surf_co2, 3=ps, 4=orb.param
     194      INTEGER :: criterion_stop  ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
    192195
    193196
    … …  
    216219!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    217220      REAL ,allocatable :: watercap_slope(:,:)                           ! Physics x Nslope: watercap per slope
     221      REAL ,allocatable :: watercap_slope_saved                          ! Value saved at the previous time step
    218222      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
    219223      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
    … …  
    229233      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
    230234      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    231       REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial co2 ice
     235      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
    232236      REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
    233237      REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
    … …  
    257261     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
    258262     REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
    259      REAL,ALLOCATABLE  :: alph_locslope(:,:)  ! Physic x Soil: Intermediate to compute Tsoil [1]
    260      REAL,ALLOCATABLE  :: beta_locslope(:,:)  ! Physic x Soil : Intermediate tocompute Tsoil [K]
    261263     REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
    262264     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
    … …  
    278280    INTEGER ::  year_iter                ! number of iteration
    279281    INTEGER :: year_iter_max             ! maximum number of iterations before stopping
    280     REAL :: timestep                     ! timestep [s]
     282    REAL :: timestep                     ! timestep [s]
     283    REAL :: watercap_sum                 ! total mass of water cap [kg/m^2]
     284
    281285    REAL WC_sum
    282286
    … …  
    424428       enddo
    425429     enddo
    426 
     430   
    427431     call surfini(ngrid,qsurf)
     432     call surfini(ngrid,qsurf)
    428433
    429434#else
    430          call phys_state_var_init(nq)
     435        call phys_state_var_init(nq)
    431436         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
    432437         call initracer(ngrid,nq)
    … …  
    615620      call end_comsoil_h_PEM
    616621      call ini_comsoil_h_PEM(ngrid,nslope)
    617 
     622      call end_adsorption_h_PEM
     623      call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    618624      allocate(ice_depth(ngrid,nslope))
    619625      ice_depth(:,:) = 0.
    … …  
    800806      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    801807      watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
    802       co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded)
    803 
    804       do ig=1,ngrid
    805          do islope=1,nslope
    806            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)
    807          enddo
     808      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
     809
     810    do ig = 1,ngrid
     811      do islope = 1,nslope
     812          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    808813      enddo
    809 
    810     if(soil_pem) then
     814    enddo
     815
     816
     817
     818    if(adsorption_pem) then
    811819     totmassco2_adsorbded = 0.
    812820     totmassh2o_adsorbded = 0.
    … …  
    878886       print *, 'Global average pressure new time step',global_ave_press_new
    879887     
    880      if(soil_pem) then
     888     if(adsorption_pem) then
    881889      do i=1,ngrid
    882890        global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    … …  
    986994
    987995!------------------------
    988 
     996     
    989997! II_d.1 Update Tsurf
    990998
    … …  
    10631071  allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
    10641072  allocate(Tsurf_locslope(ngrid))
    1065   allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
    1066   allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
    1067 
    10681073  print *,"Updating soil temperature"
    10691074
    … …  
    10711076  do islope = 1,nslope
    10721077     TI_locslope(:,:) = TI_PEM(:,:,islope)
    1073      alph_locslope(:,:) = alph_PEM(:,:,islope)
    1074      beta_locslope(:,:) = beta_PEM(:,:,islope)
    10751078     do t = 1,timelen
    10761079       Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
    10771080       Tsurf_locslope(:) =  tsurf_phys_GCM_timeseries(:,islope,t)
    1078 
    1079        call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    1080        call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    1081 
    1082 
     1081       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
     1082       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    10831083       tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
    10841084       do ig = 1,ngrid
    10851085         do isoil = 1,nsoilmx_PEM
    1086           watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
     1086          watersoil_density_phys_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
    10871087          if(isnan(Tsoil_locslope(ig,isoil))) then
    1088              write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
    1089              stop
     1088            call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)
    10901089          endif
    10911090         enddo
    10921091       enddo
    1093      enddo
    1094      alph_PEM(:,:,islope) = alph_locslope(:,:)
    1095      beta_PEM(:,:,islope) = beta_locslope(:,:)
     1092
     1093
     1094     enddo
    10961095  enddo
    1097 
    1098  
    10991096     tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
    11001097     watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
     1098
     1099
    11011100  print *, "Update of soil temperature done"
    11021101
    … …  
    11041103  deallocate(Tsoil_locslope)
    11051104  deallocate(Tsurf_locslope)
    1106   deallocate(alph_locslope)
    1107   deallocate(beta_locslope)
    1108 
    11091105  write(*,*) "Compute ice table"
    11101106
    11111107! II_d.3 Update the ice table
    1112      call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
     1108     call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
    11131109       
    11141110      print *, "Update soil propreties"
    11151111! II_d.4 Update the soil thermal properties
    1116       call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
     1112      call update_soil(ngrid,nslope,nsoilmx,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
    11171113        ice_depth,TI_PEM)
    11181114
    11191115! II_d.5 Update the mass of the regolith adsorbded
    1120 
     1116     if(adsorption_pem) then
    11211117    call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
    11221118                             tsoil_PEM,TI_PEM,ps_phys_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    11231119                             h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    1124 
     1120 
     1121     endif
    11251122  endif !soil_pem
    11261123
    … …  
    11511148
    11521149!------------------------
    1153       call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
    1154 
    1155       call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
     1150      call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
     1151
     1152      call criterion_co2_stop(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid, &
     1153                                   initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
    11561154
    11571155      year_iter=year_iter+dt_pem
     1156
     1157     
    11581158
    11591159      print *, "Checking all the stopping criterion."
    11601160      if (STOPPING_water) then
    1161         print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
     1161        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
     1162        criterion_stop=1
    11621163        criterion_stop=1
    11631164      endif
    … …  
    11831184      endif
    11841185
    1185 
    1186 
    1187       if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
     1186     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
    11881187        exit
    11891188      else
    11901189        print *, "We continue!"
    1191         print *, "Number of iteration of the PEM : year_iter=", year_iter 
     1190        print *, "Number of iteration of the PEM : year_iter=", year_iter
    11921191      endif
     1192
     1193
    11931194
    11941195      global_ave_press_old=global_ave_press_new
    … …  
    12211222#endif
    12221223      ENDDO ! of DO ig=1,ngrid
     1224
     1225
     1226! H2O ice
     1227
     1228
     1229
     1230     DO ig=1,ngrid
     1231       if(watercaptag(ig)) then
     1232       
     1233         watercap_sum=0.
     1234         DO islope=1,nslope
     1235           watercap_slope_saved = watercap_slope(ig,islope)
     1236           if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then
     1237             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     1238           else
     1239             watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     1240             qsurf_slope(ig,igcm_h2o_ice,islope)=0.
     1241           endif
     1242           watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     1243!           watercap_slope(ig,islope)=0.
     1244         enddo
     1245           water_reservoir(ig)=water_reservoir(ig)+watercap_sum
     1246       endif
     1247     enddo
     1248
     1249
     1250      DO ig = 1,ngrid
     1251        qsurf(ig,igcm_h2o_ice) = 0.
     1252        DO islope = 1,nslope
     1253          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
     1254                                   * subslope_dist(ig,islope) /          &
     1255                                  cos(pi*def_slope_mean(islope)/180.)
     1256        ENDDO
     1257      ENDDO ! of DO ig=1,ngrid
     1258
     1259     DO ig=1,ngrid
     1260!       DO islope=1,nslope
     1261         if(qsurf(ig,igcm_h2o_ice).GT.500) then
     1262            watercaptag(ig)=.true.
     1263       DO islope=1,nslope
     1264            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
     1265            water_reservoir(ig)=water_reservoir(ig)+250
     1266       ENDDO
     1267         endif
     1268!       enddo
     1269     enddo
     1270
     1271     DO ig=1,ngrid
     1272         if(water_reservoir(ig).LT. 10) then
     1273            watercaptag(ig)=.false.
     1274            qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
     1275       DO islope=1,nslope
     1276            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
     1277       ENDDO
     1278         endif
     1279
     1280     enddo
     1281
     1282      DO ig = 1,ngrid
     1283        watercap(ig) = 0.
     1284        DO islope = 1,nslope
     1285          watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
     1286                                   * subslope_dist(ig,islope) /          &
     1287                                  cos(pi*def_slope_mean(islope)/180.)
     1288        ENDDO
     1289      ENDDO ! of DO ig=1,ngrid
     1290
     1291
     1292
     1293
     1294
     1295     
     1296
    12231297
    12241298! III_a.2  Tsoil update (for startfi)
    … …  
    14601534
    14611535!------------------------
    1462 
    1463          call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
     1536        call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
     1537                         float(day_ini),0.,nslope,def_slope,subslope_dist)
     1538
     1539
     1540        call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    14641541                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    1465 
    1466      call info_run_PEM(year_iter, criterion_stop)
     1542        call info_run_PEM(year_iter, criterion_stop)
    14671543
    14681544      print *, "restartfi_PEM.nc has been written"

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

    r2885 r2888  
    11   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    3                       global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys, m_h2o_regolith_phys,deltam_h2o_regolith_phys)
     3                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
     4                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    45   
    56   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    6    use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM, water_reservoir
     7   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM,water_reservoir_nom
    78   use comsoil_h, only:  volcapa,inertiedat
    8    use adsorption_mod, only : regolith_adsorption
     9   use adsorption_mod, only : regolith_adsorption,adsorption_pem
    910   USE temps_mod_evol, ONLY: year_PEM
     11   USE ice_table_mod, only: computeice_table_equilibrium
    1012#ifndef CPP_STD   
    1113   use surfdat_h, only: watercaptag
    1214#endif
    13 
    1415
    1516  implicit none
    … …  
    4546  real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of h2o adsorbed [kg/m^2]
    4647  real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     48  real,intent(inout) :: water_reservoir(ngrid)                        ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
    4749! local
    4850   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
    … …  
    7072!!! Purpose: read start_pem. Need a specific iostart_PEM
    7173!!!
    72 !!! Order: 1. Previous year of the PEM run
    73 !!!        2. Thermal Inertia
    74 !!!        3. Soil Temperature
    75 !!!        4. Ice table
    76 !!!        5. Mass of CO2 adsorbed
     74!!! Order: 0. Previous year of the PEM run
     75!!!        1. Thermal Inertia
     76!!!        2. Soil Temperature
     77!!!        3. Ice table
     78!!!        4. Mass of CO2 & H2O adsorbed
     79!!!        5. Water reservoir
    7780!!!
    7881!!! /!\ This order must be respected !
    … …  
    8689inquire(file=filename,exist =  startpem_file)
    8790
    88 write(*,*)'Start PEM=',startpem_file
    89 
    90 
     91write(*,*)'Is start PEM?',startpem_file
     92
     93startpem_file = .false.
    9194!1. Run
    9295
    … …  
    213216            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
    214217      enddo
    215     enddo
    216 
     218          do iloop = nsoil_GCM+1,nsoil_PEM
     219             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
     220          enddo
     221      enddo
     222!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     223!     do iloop = 1,2000000
     224!      write(*,*) 'iloop,islope',islope,iloop
     225!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     226!     enddo
     227   
    217228   else
    218229! predictor corrector: restart from year 1 of the GCM and build the evolution of
    … …  
    220231
    221232    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    222     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
    223     call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
     233    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
     234    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    224235    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
    225     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
     236    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
    226237
    227238
    … …  
    243254
    244255   
    245 ENDDO
     256ENDDO ! islope
    246257
    247258print *,'PEMETAT0: SOIL TEMP  DONE'
    … …  
    253264
    254265
    255      call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
    256      call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
     266     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
     267     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
    257268
    258269
    … …  
    262273
    263274!4. CO2 & H2O Adsorption
     275 if(adsorption_pem) then
    264276  DO islope=1,nslope
    265277   write(num,fmt='(i2.2)') islope
    … …  
    290302    endif
    291303print *,'PEMETAT0: CO2 & H2O adsorption done  '
    292 
     304    endif
    293305endif ! soil_pem
     306
     307
     308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     309
     310!. 5 water reservoir
     311
     312#ifndef CPP_STD   
     313   call get_field("water_reservoir",water_reservoir,found)
     314   if(.not.found) then
     315      write(*,*)'Pemetat0: failed loading <water_reservoir>'
     316      write(*,*)'will reconstruct the values from watercaptag'
     317      do ig=1,ngrid
     318        if(watercaptag(ig)) then
     319           water_reservoir=water_reservoir_nom
     320        else
     321           water_reservoir=0.
     322        endif
     323      enddo
     324   endif
     325#endif
     326
    294327
    295328call close_startphy
    … …  
    396429            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
    397430      enddo
    398      enddo
    399    
     431     
     432           do iloop = nsoil_GCM+1,nsoil_PEM
     433             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
     434          enddo
     435       enddo
     436!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     437!     do iloop = 1,2000000
     438!      write(*,*) 'islope,iloop',islope,iloop
     439!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     440!     enddo
     441 
     442   
    400443    do it = 1,timelen
    401444        do isoil = nsoil_GCM+1,nsoil_PEM
    … …  
    412455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    413456!c) Ice table   
    414      call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
    415      call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
     457     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
     458     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
    416459
    417460
    … …  
    425468!d) Regolith adsorbed
    426469
    427 
     470 if(adsorption_pem) then
    428471    m_co2_regolith_phys(:,:,:) = 0.
    429472    m_h2o_regolith_phys(:,:,:) = 0.
    … …  
    434477    deltam_co2_regolith_phys(:) = 0.
    435478    deltam_h2o_regolith_phys(:) = 0. 
     479 endif
     480
     481print *,'PEMETAT0: CO2 adsorption done  '
     482
     483endif !soil_pem
     484
     485endif ! of if (startphy_file)
     486
     487!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    436488   
    437489
    438 print *,'PEMETAT0: CO2 adsorption done  '
    439 
    440 endif !soil_pem
    441 
    442 endif ! of if (startphy_file)
    443 
    444 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    445    
    446      if(soil_pem) then
    447 
    448    
    449 
     490
     491!. e) water reservoir
     492
     493#ifndef CPP_STD   
     494
     495      do ig=1,ngrid
     496        if(watercaptag(ig)) then
     497           water_reservoir=water_reservoir_nom
     498        else
     499           water_reservoir=0.
     500        endif
     501      enddo
     502
     503#endif
     504
     505if(soil_pem) then
    450506!! Sanity check
    451507    DO ig = 1,ngrid
    … …  
    453509        DO iloop = 1,nsoil_PEM
    454510           if(isnan(tsoil_PEM(ig,iloop,islope))) then
    455          write(*,*) "failed nan construction", ig, iloop, islope
    456            stop
     511         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1)
    457512           endif
    458513        ENDDO
    459514      ENDDO
    460515     ENDDO
    461 
    462      endif!soil_pem
     516endif!soil_pem
    463517
    464518

  • trunk/LMDZ.COMMON/libf/evolution/pemredem.F90

    r2885 r2888  
    6464  call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
    6565
    66   call put_var("FluxGeo","Geothermal Flux (W/m^2)",fluxgeo)
    6766
    6867  ! Close file
    … …  
    7473                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
    7574                    m_co2_regolith,m_h2o_regolith,water_reservoir)
     75
     76
    7677  ! write time-dependent variable to restart file
    7778  use iostart_PEM, only : open_restartphy, close_restartphy, &
    … …  
    9697  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
    9798  real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
    98   real,intent(IN) :: water_reservoir(ngrid)
     99  real,intent(in) :: water_reservoir(ngrid)
    99100  integer :: islope
    100101  CHARACTER*2 :: num 
    … …  
    116117  ! set time counter in file
    117118  call put_var("Time","Year of simulation",year_tot)
    118 
     119  call put_field('water_reservoir','water_reservoir', &
     120           water_reservoir,year_tot)
    119121  call put_field("water_reservoir","water_reservoir", &
    120122                 water_reservoir,year_tot)

  • trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90

    r2885 r2888  
    6969  REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                               ! co2 ice, mesh averaged, of the concatenated file [kg/m^2]
    7070  REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! h2o ice per slope of the concatenated file [kg/m^2]
    71   REAL, ALLOCATABLE :: watercap_slope(:,:,:,:)
    72 
     71  REAL, ALLOCATABLE ::  watercap_slope(:,:,:,:)
    7372!-----------------------------------------------------------------------
    7473  modname="read_data_gcm"
    … …  
    8382      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
    8483      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
     84      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
    8585
    8686  print *, "Opening ", fichnom, "..."
    … …  
    135135
    136136     print *, "Downloading data for h2o_ice_s_slope done"
     137
    137138     print *, "Downloading data for watercap_slope ..."
    138 
    139 DO islope=1,nslope
    140   write(num,fmt='(i2.2)') islope
    141   call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
    142 ENDDO
    143 
     139DO islope=1,nslope
     140       write(num,fmt='(i2.2)') islope
     141       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))   
     142ENDDO           
    144143     print *, "Downloading data for watercap_slope done"
    145      print *, "Downloading data for tsurf_slope ..."
     144   
     145 print *, "Downloading data for tsurf_slope ..."
    146146
    147147DO islope=1,nslope

  • trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90

    r2855 r2888  
    11      subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
    2                therm_i,                          &
    3                timestep,tsurf,tsoil,alph_PEM,beta_PEM)
     2               therm_i,timestep,tsurf,tsoil)
    43
    54
    65      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
    76                          mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
    8                           coefd_PEM, mu_PEM,fluxgeo
     7                          coefd_PEM, mu_PEM,alph_PEM,beta_PEM,fluxgeo
    98      use comsoil_h,only: volcapa
    109      implicit none
    … …  
    3433! outputs:
    3534      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
    36       real,intent(inout) :: alph_PEM(ngrid,nsoil-1) ! intermediate for computations [1]
    37       real,intent(inout) :: beta_PEM(ngrid,nsoil-1) ! intermediate for computations [K]
    38      
    3935! local variables:
    4036      integer ig,ik   

  • trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90

    r2835 r2888  
    55  INTEGER   year_bp_ini     !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
    66  INTEGER   dt_pem          !     dt_pem  : in years, the time step used by the PEM   
    7   REAL      alpha_criterion !     alpha_criterion : percentage of change before stopping the PEM
     7  REAL      ice_criterion   !     ice_criterion : percentage of change of ice before stopping the PEM
     8  REAL      ps_criterion    !     ice_criterion : percentage of change of ice before stopping the PEM
    89  INTEGER   year_PEM        !     year written in startfiPEM.nc
    910  INTEGER   Max_iter_pem    !     Maximal number of iteration when converging to a steady state, read in evol.def

  • trunk/LMDZ.COMMON/libf/evolution/update_soil.F90

    r2863 r2888  
    1    SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,&
    2                           ice_depth,TI_PEM)
     1   SUBROUTINE update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,TI_PEM)
    32#ifndef CPP_STD
    43 USE comsoil_h, only:  inertiedat, volcapa
    … …  
    87 implicit none
    98! Input:
    10  INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
    11  REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
     9 INTEGER,INTENT(IN) ::  ngrid, nslope,nsoil_GCM, nsoil_PEM
    1210 REAL,INTENT(IN) :: p_avg_new
    13  REAL,INTENT(IN) :: co2ice(ngrid,nslope)
     11 REAL,INTENT(IN) :: tendencies_waterice(ngrid,nslope)
    1412 REAL,INTENT(IN) :: waterice(ngrid,nslope)
    1513 REAL,INTENT(in) :: ice_depth(ngrid,nslope)
    16 
    1714! Outputs:
    1815
    … …  
    2118! Constants:
    2219
    23  REAL ::  alpha = 0.2
    24  REAL ::  beta = 1.08e7
    25  REAL ::  To = 273.15
    26  REAL ::  R =  8.314
    27  REAL ::  L =  51058.
    2820 REAL ::  inertie_thresold = 800. ! look for ice
    2921 REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
    3022 REAL ::  ice_inertia = 1200  ! Inertia of ice
    31  REAL ::  P610 = 610.0
    32  REAL ::  m_h2o = 18.01528E-3
    33  REAL ::  m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
    34  REAL ::  m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
    35  REAL ::  A,B,mmean             
     23 REAL ::  P610 = 610.0       ! current average pressure on Mars [Pa]
     24 REAL :: TI_breccia = 750.   ! THermal inertia of Breccia following Wood 2009 [SI]
     25 REAL :: TI_bedrock = 2300.  ! Thermal inertia of Bedrock following Wood 2009 [SI]
    3626
    3727! Local variables:
    … …  
    3929 INTEGER :: ig,islope,iloop,iref,k
    4030 REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith
    41  REAL :: d(ngrid,nsoil_PEM,nslope)
    42  REAL :: Total_surface
    43  INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
    44  INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
    45 
    46 ! 0. Initialisation
    47 
    48 !  do ig = 1,ngrid
    49 !    do islope = 1,nslope
    50 !     if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then
    51 !        ispermanent_h2oglaciers(ig,islope) = 1
    52 !     else
    53 !        ispermanent_h2oglaciers(ig,islope) = 0
    54 !     endif
    55 !
    56 !     if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then
    57 !        ispermanent_co2glaciers(ig,islope) = 1
    58 !     else
    59 !        ispermanent_co2glaciers(ig,islope) = 0
    60 !     endif
    61 !    enddo
    62 !  enddo
    63 
    64 
    65 
     31 REAL :: delta
     32 REAL :: TI_breccia_new
    6633!  1.Ice TI feedback
    6734
    … …  
    7239! 2. Modification of the regolith thermal inertia.
    7340
    74 !  do ig=1,ngrid
    75 !    do islope=1,nslope
    76 !     do iloop =1,n_1km
    77 !       d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
    78 !   if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
    79 !          TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))
    80 !
    81 !   endif
    82 !     enddo
    83 !   enddo
    84 !enddo
    8541
    86 !  3. Build new TI for the PEM
     42
     43
     44 do islope = 1,nslope
     45   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
     46   do ig = 1,ngrid
     47
     48      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
     49              regolith_inertia(ig,islope) = inertie_averaged
     50       endif
     51        write(*,*) 'ig,islope',ig,islope,inertie_thresold,TI_PEM(ig,1,islope)
     52       if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
     53!          regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
     54       endif
     55       TI_breccia_new = TI_breccia !*(p_avg_new/P610)**0.3
     56   enddo
     57 enddo
     58
     59
     60! 3. Build new Thermal Inertia
     61
     62do  islope=1,nslope
     63   do ig = 1,ngrid
     64     do iloop = 1,nsoil_GCM
     65         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
     66     enddo
     67     if(regolith_inertia(ig,islope).lt.TI_breccia_new) then
     68!!! transition
     69             delta = 50.
     70             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
     71            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
     72                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia_new**2))))
     73            do iloop=nsoil_GCM+2,n_1km
     74              TI_PEM(ig,iloop,islope) = TI_breccia_new
     75            enddo     
     76      else ! we keep the high ti values
     77            do iloop=nsoil_GCM+1,n_1km
     78              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
     79            enddo
     80       endif ! TI PEM and breccia comparison
     81!! transition
     82       delta = 1000.
     83       TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
     84            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
     85            ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
     86       do iloop=n_1km+2,nsoil_PEM
     87            TI_PEM(ig,iloop,islope) = TI_bedrock
     88       enddo   
     89      enddo ! ig
     90ENDDO ! islope
     91
     92!  4. Build new TI in case of ice table
    8793! a) For the regolith
    8894  do ig=1,ngrid
Note: See TracChangeset for help on using the changeset viewer.