Ignore:
Timestamp:
Nov 14, 2025, 5:26:48 PM (4 weeks ago)
Author:
jbclement
Message:

PEM:
Clearing all the tags and code related to the Generic PCM.
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3571 r3961  
    7171                            tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
    7272
    73 END SUBROUTINE
     73END SUBROUTINE regolith_adsorption
    7474
    7575!=======================================================================
    … …  
    8989use vertical_layers_mod,   only: ap, bp
    9090use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
    91 
    92 #ifndef CPP_STD
    93     use comcstfi_h,   only: pi
    94 #else
    95     use comcstfi_mod, only: pi
    96 #endif
     91use comcstfi_h,            only: pi
    9792
    9893implicit none
    … …  
    147142ispermanent_co2glaciers = .false.
    148143
    149 #ifndef CPP_STD
    150     ! 0.1 Look at perennial ice
     144! 0.1 Look at perennial ice
     145do ig = 1,ngrid
     146    do islope = 1,nslope
     147        if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
     148        if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     149    enddo
     150enddo
     151
     152! 0.2 Compute the partial pressure of vapor
     153! a. the molecular mass into the column
     154do ig = 1,ngrid
     155    mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B)
     156enddo
     157
     158! b. pressure level
     159do it = 1,timelen
    151160    do ig = 1,ngrid
    152         do islope = 1,nslope
    153             if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
    154             if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     161        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
     162    enddo
     163enddo
     164
     165! c. Vapor pressure
     166pvapor = mass_mean/m_h2o*q_h2o*zplev_mean
     167pvapor_avg = sum(pvapor,2)/timelen
     168deallocate(pvapor,zplev_mean,mass_mean)
     169
     170! 1. we compute the mass of H2O adsorded in each layer of the meshes
     171do ig = 1,ngrid
     172    do islope = 1,nslope
     173        do iloop = 1,index_breccia
     174            K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
     175            if (TI_PEM(ig,iloop,islope) < inertie_thresold) then
     176                theta_h2o_adsorbed(ig,iloop,islope) = (K*pvapor_avg(ig)/(1. + K*pvapor_avg(ig)))**mu
     177            else
     178                p_sat = exp(beta_clap_h2o/tsoil_PEM(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
     179                theta_h2o_adsorbed(ig,iloop,islope) = (K*p_sat/(1. + K*p_sat))**mu
     180            endif
     181            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbed(ig,iloop,islope)*m_theta*rho_regolith
    155182        enddo
    156183    enddo
    157 
    158     ! 0.2 Compute the partial pressure of vapor
    159     ! a. the molecular mass into the column
    160     do ig = 1,ngrid
    161         mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B)
    162     enddo
    163 
    164     ! b. pressure level
    165     do it = 1,timelen
    166         do ig = 1,ngrid
    167             zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
     184enddo
     185
     186! 2. Check the exchange between the atmosphere and the regolith
     187do ig = 1,ngrid
     188    delta_mreg(ig) = 0.
     189    do islope = 1,nslope
     190        deltam_reg_slope(ig,islope) = 0.
     191        do iloop = 1,index_breccia
     192            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
     193                if (iloop == 1) then
     194                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
     195                else
     196                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
     197                endif
     198            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
     199                deltam_reg_complete(ig,iloop,islope) = 0.
     200            endif
     201            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    168202        enddo
    169     enddo
    170 
    171     ! c. Vapor pressure
    172     pvapor = mass_mean/m_h2o*q_h2o*zplev_mean
    173     pvapor_avg = sum(pvapor,2)/timelen
    174 #endif
    175 deallocate(pvapor,zplev_mean,mass_mean)
    176 
    177 #ifndef CPP_STD
    178     ! 1. we compute the mass of H2O adsorded in each layer of the meshes
    179     do ig = 1,ngrid
    180         do islope = 1,nslope
    181             do iloop = 1,index_breccia
    182                 K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
    183                 if (TI_PEM(ig,iloop,islope) < inertie_thresold) then
    184                     theta_h2o_adsorbed(ig,iloop,islope) = (K*pvapor_avg(ig)/(1. + K*pvapor_avg(ig)))**mu
    185                 else
    186                     p_sat = exp(beta_clap_h2o/tsoil_PEM(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
    187                     theta_h2o_adsorbed(ig,iloop,islope) = (K*p_sat/(1. + K*p_sat))**mu
    188                 endif
    189                 dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbed(ig,iloop,islope)*m_theta*rho_regolith
    190             enddo
    191         enddo
    192     enddo
    193 
    194     ! 2. Check the exchange between the atmosphere and the regolith
    195     do ig = 1,ngrid
    196         delta_mreg(ig) = 0.
    197         do islope = 1,nslope
    198             deltam_reg_slope(ig,islope) = 0.
    199             do iloop = 1,index_breccia
    200                 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
    201                     if (iloop == 1) then
    202                         deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
    203                     else
    204                         deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
    205                     endif
    206                 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
    207                     deltam_reg_complete(ig,iloop,islope) = 0.
    208                 endif
    209                 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    210             enddo
    211             delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    212         enddo
    213     enddo
    214     m_h2o_completesoil = dm_h2o_regolith_slope
    215 #endif
    216 END SUBROUTINE
     203        delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     204    enddo
     205enddo
     206m_h2o_completesoil = dm_h2o_regolith_slope
     207
     208END SUBROUTINE regolith_h2oadsorption
    217209
    218210!=======================================================================
    … …  
    230222use vertical_layers_mod,   only: ap, bp
    231223use constants_marspem_mod, only: m_co2, m_noco2, rho_regolith
    232 
    233 #ifndef CPP_STD
    234     use comcstfi_h,   only: pi
    235 #else
    236     use comcstfi_mod, only: pi
    237 #endif
     224use comcstfi_h,            only: pi
    238225
    239226implicit none
    … …  
    288275ispermanent_co2glaciers = .false.
    289276
    290 #ifndef CPP_STD
    291     ! 0.1 Look at perennial ice
     277! 0.1 Look at perennial ice
     278do ig = 1,ngrid
     279    do islope = 1,nslope
     280        if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
     281        if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     282    enddo
     283enddo
     284
     285! 0.2  Compute the partial pressure of CO2
     286! a. the molecular mass into the column
     287do ig = 1,ngrid
     288    mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B)
     289enddo
     290
     291! b. pressure level
     292do it = 1,timelen
    292293    do ig = 1,ngrid
    293         do islope = 1,nslope
    294             if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
    295             if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     294        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
     295    enddo
     296enddo
     297
     298! c. Vapor pressure
     299pco2 = mass_mean/m_co2*q_co2*zplev_mean
     300pco2_avg(:) = sum(pco2(:,:),2)/timelen
     301
     302deallocate(zplev_mean,mass_mean,pco2)
     303
     304! 1. Compute the fraction of the pores occupied by H2O
     305call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     306                            theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
     307
     308! 2. we compute the mass of co2 adsorded in each layer of the meshes
     309do ig = 1,ngrid
     310    do islope = 1,nslope
     311        do iloop = 1,index_breccia
     312            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
     313                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
     314                                                         (alpha*pco2_avg(ig) + sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
     315            else
     316                if (abs(m_co2_completesoil(ig,iloop,islope)) < 1.e-10) then !!! we are at first call
     317                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
     318                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
     319                else ! no change: permanent ice stick the atoms of CO2
     320                    dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
     321                endif
     322            endif
    296323        enddo
    297324    enddo
    298 
    299     ! 0.2  Compute the partial pressure of CO2
    300     ! a. the molecular mass into the column
    301     do ig = 1,ngrid
    302         mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B)
    303     enddo
    304 
    305     ! b. pressure level
    306     do it = 1,timelen
    307         do ig = 1,ngrid
    308             zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
     325enddo
     326
     327! 3. Check the exchange between the atmosphere and the regolith
     328do ig = 1,ngrid
     329    delta_mreg(ig) = 0.
     330    do islope = 1,nslope
     331        deltam_reg_slope(ig,islope) = 0.
     332        do iloop = 1,index_breccia
     333            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
     334                if (iloop == 1) then
     335                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
     336                else
     337                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
     338                endif
     339            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
     340                deltam_reg_complete(ig,iloop,islope) = 0.
     341            endif
     342            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    309343        enddo
    310     enddo
    311 
    312     ! c. Vapor pressure
    313     pco2 = mass_mean/m_co2*q_co2*zplev_mean
    314     pco2_avg(:) = sum(pco2(:,:),2)/timelen
    315 
    316     deallocate(zplev_mean,mass_mean,pco2)
    317 
    318     ! 1. Compute the fraction of the pores occupied by H2O
    319     call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    320                                 theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
    321 
    322     ! 2. we compute the mass of co2 adsorded in each layer of the meshes
    323     do ig = 1,ngrid
    324         do islope = 1,nslope
    325             do iloop = 1,index_breccia
    326                 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
    327                     dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
    328                                                              (alpha*pco2_avg(ig) + sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
    329                 else
    330                     if (abs(m_co2_completesoil(ig,iloop,islope)) < 1.e-10) then !!! we are at first call
    331                         dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
    332                                                                  /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
    333                     else ! no change: permanent ice stick the atoms of CO2
    334                         dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
    335                     endif
    336                 endif
    337             enddo
    338         enddo
    339     enddo
    340 
    341     ! 3. Check the exchange between the atmosphere and the regolith
    342     do ig = 1,ngrid
    343         delta_mreg(ig) = 0.
    344         do islope = 1,nslope
    345             deltam_reg_slope(ig,islope) = 0.
    346             do iloop = 1,index_breccia
    347                 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
    348                     if (iloop == 1) then
    349                         deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
    350                     else
    351                         deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
    352                     endif
    353                 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
    354                     deltam_reg_complete(ig,iloop,islope) = 0.
    355                 endif
    356                 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    357             enddo
    358             delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    359         enddo
    360     enddo
    361     m_co2_completesoil = dm_co2_regolith_slope
    362 #endif
     344        delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     345    enddo
     346enddo
     347m_co2_completesoil = dm_co2_regolith_slope
    363348
    364349END SUBROUTINE regolith_co2adsorption
Note: See TracChangeset for help on using the changeset viewer.