Changeset 3320 for trunk


Ignore:
Timestamp:
Apr 30, 2024, 12:08:09 PM (7 months ago)
Author:
jbclement
Message:

PEM:
Small correction of a bug in the reading of "run_PEM.def" + some cleanings.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r3317 r3320  
    286286== 25/04/2024 == JBC
    287287The "start" and "startfi" file names are renamed to match those from the Mars PCM. This is necessary to initialize correctly the 3D PEM with "phys_state_var_init_mod.F90".
     288
     289== 26/04/2024 == JBC
     290- Small corrections to make the PEM work in 3D.
     291- Addition of the flag "layering_algo" (Default = .false.) defined in the "run_PEM.def" to choose to use the layering algorithm.
     292
     293== 30/04/2024 == JBC
     294Small correction of a bug in the reading of "run_PEM.def" + some cleanings.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3319 r3320  
    9292# co2ice_flow=.true.
    9393
    94 !#---------- Layering parameters ----------#
     94#---------- Layering parameters ----------#
    9595# Do you want to run with the layering algorithm? Default = .false.
    9696# layering=.true.

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

    r3264 r3320  
    1 module ice_table_mod
    2 
    3   implicit none
     1MODULE ice_table_mod
     2
     3implicit none
    44
    55!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    1010!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1111
    12   LOGICAL,SAVE :: icetable_equilibrium                    ! Boolean to say if the PEM needs to recompute the icetable depth when at  equilibrium
    13   LOGICAL,SAVE :: icetable_dynamic                        ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method)
    14   real,save,allocatable  :: porefillingice_depth(:,:)    ! ngrid x nslope: Depth of the ice table [m]
    15   real,save,allocatable  :: porefillingice_thickness(:,:)    ! ngrid x nslope: Thickness of the ice table [m]
    16 
     12logical, save                           :: icetable_equilibrium     ! Boolean to say if the PEM needs to recompute the icetable depth when at  equilibrium
     13logical, save                           :: icetable_dynamic         ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method)
     14real, allocatable, dimension(:,:), save :: porefillingice_depth     ! ngrid x nslope: Depth of the ice table [m]
     15real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m]
     16
     17!----------------------------------------
    1718  contains
    18 
    19 
    20   subroutine ini_ice_table_porefilling(ngrid,nslope)
    21  
    22     implicit none
    23     integer,intent(in) :: ngrid ! number of atmospheric columns
    24     integer,intent(in) :: nslope ! number of slope within a mesh
    25 
    26     allocate(porefillingice_depth(ngrid,nslope))
    27     allocate(porefillingice_thickness(ngrid,nslope))
    28 
    29   end subroutine ini_ice_table_porefilling
    30 
    31   subroutine end_ice_table_porefilling
    32 
    33     implicit none
    34     if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
    35     if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
    36 
    37   end subroutine end_ice_table_porefilling
    38 
    39 !!! --------------------------------------
    40 
    41    SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
     19!----------------------------------------
     20SUBROUTINE ini_ice_table_porefilling(ngrid,nslope)
     21
     22implicit none
     23
     24integer, intent(in) :: ngrid  ! number of atmospheric columns
     25integer, intent(in) :: nslope ! number of slope within a mesh
     26
     27allocate(porefillingice_depth(ngrid,nslope))
     28allocate(porefillingice_thickness(ngrid,nslope))
     29
     30END SUBROUTINE ini_ice_table_porefilling
     31
     32!----------------------------------------
     33SUBROUTINE end_ice_table_porefilling
     34
     35implicit none
     36
     37if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
     38if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
     39
     40 END SUBROUTINE end_ice_table_porefilling
     41
     42!----------------------------------------
     43SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
    4244!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4345!!!
    … …  
    4547!!! density at the surface and at depth.
    4648!!! Computations are made following the methods in Schorgofer et al., 2005
    47 !!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
    48 !!! 
     49!!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere
     50!!!
    4951!!! Author: LL
    5052!!!
    5153!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    52     use math_mod,only: findroot
    53     USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM                  ! Depth of the vertical grid
    54     USE soil_thermalproperties_mod, only: ice_thermal_properties
    55     implicit none
    56 !   inputs
    57 ! -----------
    58     integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    59     logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
    60     real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
    61     real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
    62     real,intent(in) :: regolith_inertia(ngrid,nslope)              ! Thermal inertia of the regolith layer [SI]
     54use math_mod,                   only: findroot
     55use comsoil_h_PEM,              only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
     56use soil_thermalproperties_mod, only: ice_thermal_properties
     57
     58implicit none
     59
     60!   Inputs
     61! --------
     62integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM
     63logical, dimension(ngrid),               intent(in) :: watercaptag              ! Boolean to check the presence of a perennial glacier
     64real, dimension(ngrid,nslope),           intent(in) :: rhowatersurf_ave         ! Water density at the surface, yearly averaged [kg/m^3]
     65real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave         ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
     66real, dimension(ngrid,nslope),           intent(in) :: regolith_inertia         ! Thermal inertia of the regolith layer [SI]
    6367!   Ouputs
    64 ! -----------
    65     real,intent(out) :: ice_table_beg(ngrid,nslope)               ! ice table depth [m]
    66     real,intent(out) :: ice_table_thickness(ngrid,nslope)        ! ice table thickness [m]
    67 !   Local
    68 ! -----------
    69     integer ig, islope,isoil,isoilend                              ! loop variables
    70     real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    71     real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
    72     real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
    73     real :: stretch                                                ! stretch factor to improve the convergence of the ice table
    74     real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]   
     68! --------
     69real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
     70real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
     71!   Locals
     72! --------
     73integer ig, islope,isoil,isoilend                              ! loop variables
     74real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
     75real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
     76real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
     77real :: stretch                                                ! stretch factor to improve the convergence of the ice table
     78real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]
    7579!   Code
    76 ! -----------
    77 
     80! ------
    7881     previous_icetable_depth(:,:) = ice_table_beg(:,:)
    79      do ig = 1,ngrid
    80       if(watercaptag(ig)) then
    81          ice_table_beg(ig,:) = 0.
    82          ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    83       else
     82do ig = 1,ngrid
     83    if (watercaptag(ig)) then
     84        ice_table_beg(ig,:) = 0.
     85        ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     86    else
    8487        do islope = 1,nslope
    85            ice_table_beg(ig,islope) = -1.
    86            ice_table_thickness(ig,islope) = 0.
    87          do isoil = 1,nsoil_PEM
    88            diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
    89          enddo
    90          if(diff_rho(1) > 0) then                                  ! ice is at the surface
    91            ice_table_beg(ig,islope) = 0.
    92                do isoilend = 2,nsoil_PEM-1
    93                  if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
    94                   call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
    95                   ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
    96                   exit
    97                  endif
    98                enddo
    99          else
    100            do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
    101              if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
    102                call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope))
    103                ! Now let's find the end of the ice table
    104                ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    105                do isoilend = isoil+1,nsoil_PEM-1
    106                  if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
    107                   call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
    108                   ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
    109                   exit
    110                  endif
    111                enddo
    112              endif !diff_rho(z) <0 & diff_rho(z+1) > 0
    113            enddo
    114           endif ! diff_rho(1) > 0
     88            ice_table_beg(ig,islope) = -1.
     89            ice_table_thickness(ig,islope) = 0.
     90            do isoil = 1,nsoil_PEM
     91                diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
     92            enddo
     93            if (diff_rho(1) > 0) then ! ice is at the surface
     94                ice_table_beg(ig,islope) = 0.
     95                do isoilend = 2,nsoil_PEM-1
     96                    if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend+1) < 0.) then
     97                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
     98                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
     99                        exit
     100                    endif
     101                enddo
     102            else
     103                do isoil = 1,nsoil_PEM - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
     104                    if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then
     105                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope))
     106                        ! Now let's find the end of the ice table
     107                        ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     108                        do isoilend = isoil+1,nsoil_PEM-1
     109                            if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
     110                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
     111                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
     112                                exit
     113                            endif
     114                        enddo
     115                    endif !diff_rho(z) <0 & diff_rho(z+1) > 0
     116                enddo
     117            endif ! diff_rho(1) > 0
    115118        enddo
    116        endif ! watercaptag
    117       enddo
     119    endif ! watercaptag
     120enddo
    118121
    119122! Small trick to speed up the convergence, Oded's idea.
    120       do islope = 1,nslope
    121         do ig = 1,ngrid
    122          if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then
     123do islope = 1,nslope
     124    do ig = 1,ngrid
     125        if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
    123126            call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
    124127            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
    125 
    126128            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
    127             previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
    128             ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
    129          endif
    130         enddo
    131       enddo
    132 
    133       RETURN
    134       END
    135 
    136 !!! --------------------------------------
    137 
    138 
    139    SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
    140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    141 !!!
    142 !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
    143 !!!
     129                                             previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
     130            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
     131        endif
     132    enddo
     133enddo
     134
     135END SUBROUTINE
     136
     137!----------------------------------------
     138SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
     139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     140!!!
     141!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
     142!!!
    144143!!! Author: LL
    145144!!!
    146145!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    147     use comsoil_h_PEM, only: mlayer_PEM
    148     use comslope_mod, only: subslope_dist,def_slope_mean
    149     use constants_marspem_mod, only: porosity
    150     use vdifc_mod, only: compute_Tice
     146use comsoil_h_PEM,        only: mlayer_PEM
     147use comslope_mod,          only: subslope_dist, def_slope_mean
     148use constants_marspem_mod, only: porosity
     149use vdifc_mod,            only: compute_Tice
    151150#ifndef CPP_STD
    152151    use comcstfi_h,   only: pi
    … …  
    155154#endif
    156155
    157     implicit none
    158 !   inputs
    159 ! -----------
    160     integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    161     real,intent(in) :: former_ice_table_thickness(ngrid,nslope)    ! ice table thickness at the former iteration [m]
    162     real,intent(in) :: new_ice_table_thickness(ngrid,nslope)       ! ice table thickness at the current iteration [m]
    163     real,intent(in) :: ice_table_depth(ngrid,nslope)              ! ice table depth [m]
    164     real,intent(in) :: tsurf(ngrid,nslope)                         ! Surface temperature [K]
    165     real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope)               ! Soil temperature [K]
    166 !   outputs
    167 ! -----------
    168     real,intent(out) :: delta_m_h2o(ngrid)                         ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
    169 
    170 ! local
    171     real :: rho(ngrid,nslope)                                      ! density of water ice [kg/m^3]
    172     integer :: ig,islope,ilay,iref                                 ! loop index 
    173     real :: Tice(ngrid,nslope)                                     ! ice temperature [k]
     156implicit none
     157
     158!   Inputs
     159! --------
     160integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
     161real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
     162real, dimension(ngrid,nslope),           intent(in) :: new_ice_table_thickness    ! ice table thickness at the current iteration [m]
     163real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth            ! ice table depth [m]
     164real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
     165real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
     166!   Outputs
     167! ---------
     168real, dimension(ngrid), intent(out) :: delta_m_h2o ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
     169!   Locals
     170!---------
     171integer :: ig, islope, ilay, iref     ! loop index
     172real, dimension(ngrid,nslope) :: rho  ! density of water ice [kg/m^3]
     173real, dimension(ngrid,nslope) :: Tice ! ice temperature [k]
    174174!   Code
    175 ! -----------
    176    rho(:,:) = 0.
    177    Tice(:,:) = 0.
     175! ------
     176rho = 0.
     177Tice = 0.
    178178!1. First let's compute Tice using a linear interpolation between the mlayer level
    179    do ig = 1,ngrid
    180       do islope = 1,nslope
    181          call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope), ice_table_depth(ig,islope), Tice(ig,islope))
    182          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
    183       enddo
    184    enddo     
    185 
    186          
     179do ig = 1,ngrid
     180    do islope = 1,nslope
     181        call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
     182        rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
     183    enddo
     184enddo
     185
    187186!2. Let's compute the amount of ice that has sublimated in each subslope
    188    do ig = 1,ngrid
    189       do islope = 1,nslope
    190         delta_m_h2o(ig)  = delta_m_h2o(ig) +  porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses 
    191                            *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    192       enddo
    193    enddo
    194 
    195    return
    196 end subroutine
    197 
    198 !!! --------------------------------------
    199 
    200    SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
    201 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    202 !!!
    203 !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
    204 !!!
     187do ig = 1,ngrid
     188    do islope = 1,nslope
     189        delta_m_h2o(ig) = delta_m_h2o(ig) +  porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses
     190                          *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
     191    enddo
     192enddo
     193
     194END SUBROUTINE
     195
     196!----------------------------------------
     197
     198SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
     199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     200!!!
     201!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
     202!!!
    205203!!! Author: LL
    206204!!!
    207205!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    208 use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
    209 use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple
    210 implicit none
    211 ! inputs
     206use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM
     207use math_mod,      only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
     208
     209implicit none
     210
     211! Inputs
    212212! ------
    213 
    214 real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
    215 real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
    216 real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
    217 real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
    218 real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
    219 real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
    220 integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
     213real, dimension(nsoilmx_PEM), intent(in) :: porefill    ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
     214real, dimension(nsoilmx_PEM), intent(in) :: psat_soil   ! Soil water pressure at saturation, yearly averaged [Pa]
     215real,                         intent(in) :: psat_surf   ! surface water pressure at saturation, yearly averaged [Pa]
     216real,                         intent(in) :: pwat_surf   ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
     217real,                         intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
     218real,                         intent(in) :: B           ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
     219integer,                      intent(in) :: index_IS    ! index of the soil layer where the ice sheet begins [1]
    221220real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
    222 
    223 ! outputs
     221! Outputs
    224222! -------
    225 integer,intent(out) :: index_filling ! index where the pore filling begins [1]
    226 integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
    227 real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
    228 real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
    229 
    230 ! local
    231 
    232 real :: eta(nsoilmx_PEM)  ! constriction
    233 integer :: ilay ! index for loop
    234 real :: old_depth_filling ! depth_filling saved
    235 real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
    236 integer :: index_tmp ! for loop
    237 real :: Jdry ! flux trought the dry layer
    238 real :: Jsat ! flux trought the ice layer
    239 real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
    240 integer :: index_firstice ! first index where ice appears (i.e., f > 0)
    241 real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
    242 real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
    243 real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
    244 
    245 ! constant
    246 real :: pvap2rho = 18.e-3/8.314
    247 !
    248  
     223integer,                      intent(out) :: index_filling    ! index where the pore filling begins [1]
     224integer,                      intent(out) :: index_geothermal ! index where the ice table stops [1]
     225real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
     226real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
     227!   Locals
     228!---------
     229real, dimension(nsoilmx_PEM) :: eta                          ! constriction
     230real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
     231real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
     232real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
     233integer                      :: ilay, index_tmp              ! index for loop
     234real                         :: old_depth_filling            ! depth_filling saved
     235real                         :: Jdry                         ! flux trought the dry layer
     236real                         :: Jsat                         ! flux trought the ice layer
     237real                         :: Jdry_prevlay, Jsat_prevlay   ! same but for the previous ice layer
     238integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
     239real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
     240!   Constant
     241!-----------
     242real, parameter :: pvap2rho = 18.e-3/8.314
     243!   Code
     244!-------
    249245! 0. Compute constriction over the layer
    250246! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
    251247if (index_IS.lt.0) then
    252    index_tmp = nsoilmx_PEM
    253    do ilay = 1,nsoilmx_PEM
    254        call constriction(porefill(ilay),eta(ilay))
    255    enddo
     248    index_tmp = nsoilmx_PEM
     249    do ilay = 1,nsoilmx_PEM
     250        call constriction(porefill(ilay),eta(ilay))
     251    enddo
    256252else
    257    index_tmp = index_IS
    258    do ilay = 1,index_IS-1
    259        call constriction(porefill(ilay),eta(ilay))
    260    enddo
    261    do ilay = index_IS,nsoilmx_PEM
    262        eta(ilay) = 0.
    263    enddo
     253    index_tmp = index_IS
     254    do ilay = 1,index_IS-1
     255        call constriction(porefill(ilay),eta(ilay))
     256    enddo
     257    do ilay = index_IS,nsoilmx_PEM
     258        eta(ilay) = 0.
     259    enddo
    264260endif
    265261
    266 ! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)   
    267 
     262! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
    268263old_depth_filling = depth_filling
    269264
    270 call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 
     265call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)
    271266
    272267do ilay = 1,index_tmp
    273   Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
    274   Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
    275   if((Jdry - Jsat).le.0) then
    276      index_filling = ilay
    277      exit
    278    endif
    279 enddo
    280 
    281 if(index_filling.eq.1)  depth_filling = mlayer_PEM(1)
    282 
    283 if(index_filling.gt.1) then
    284     Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
    285     Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
    286    call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
     268    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
     269    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
     270    if (Jdry - Jsat <= 0) then
     271        index_filling = ilay
     272        exit
     273    endif
     274enddo
     275
     276if (index_filling == 1) then
     277    depth_filling = mlayer_PEM(1)
     278else if (index_filling > 1) then
     279    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1)
     280    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
     281    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling)
    287282endif
    288283
    289 
    290284! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
    291 
    292 
    293285! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
    294 
    295286index_firstice = -1
    296287do ilay = 1,nsoilmx_PEM
    297    if (porefill(ilay).le.0.) then
     288    if (porefill(ilay) <= 0.) then
    298289        index_firstice = ilay  ! first point with ice
    299290        exit
    … …  
    302293
    303294! 2.1: now we can computeCompute d_z (eta* d_z(rho))
    304 
    305 call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
    306 
    307 if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
    308      call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
     295call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
     296
     297if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
     298
     299call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
     300dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
     301
     302! 3. Ice table boundary due to geothermal heating
     303if (index_IS > 0) index_geothermal = -1
     304if (index_geothermal < 0) depth_geothermal = -1.
     305if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
     306    index_geothermal = -1
     307    do ilay = 2,nsoilmx_PEM
     308        if (dz_psat(ilay) > 0.) then  ! first point with reversed flux
     309            index_geothermal = ilay
     310            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
     311            exit
     312        endif
     313    enddo
     314else
     315    index_geothermal = -1
    309316endif
    310 
    311 call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
    312 dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
    313 
    314 ! 3. Ice table boundary due to geothermal heating
    315 
    316 if(index_IS.gt.0) index_geothermal = -1
    317 if(index_geothermal.lt.0) depth_geothermal = -1.
    318 if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
    319    index_geothermal = -1
    320    do ilay=2,nsoilmx_PEM
    321         if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
    322            index_geothermal=ilay
    323            call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
    324            exit
    325         endif
    326      enddo
    327   else
    328      index_geothermal = -1
    329   endif
    330  if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
    331      call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
    332      index_tmp = -1
    333      do ilay=index_geothermal,nsoilmx_PEM
    334         if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
    335          call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
    336         if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
    337            if (ilay>index_geothermal) then
    338 !              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
    339               call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 
    340                    dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
    341 !               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*)
    342 !                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
    343               index_tmp=ilay
     317if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
     318    call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
     319    index_tmp = -1
     320    do ilay = index_geothermal,nsoilmx_PEM
     321        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full
     322        call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
     323        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
     324            if (ilay > index_geothermal) then
     325!                write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
     326                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
     327                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
     328!                if (depth_geothermal > mlayer_PEM(ilay) .or. depth_geothermal < mlayer_PEM(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay - 1),depth_geothermal,mlayer_PEM(ilay)
     329              index_tmp = ilay
    344330           endif
    345331           ! otherwise leave depth_geothermal unchanged
    … …  
    347333        endif
    348334        massfillabove = massfillafter
    349      enddo
    350      if (index_tmp.gt.0) index_geothermal = index_tmp
    351   end if
    352   return
    353 end subroutine
    354 
    355 !!! --------------------------------------
    356 
     335    enddo
     336    if (index_tmp > 0) index_geothermal = index_tmp
     337end if
     338
     339END SUBROUTINE
     340
     341!----------------------------------------
    357342SUBROUTINE constriction(porefill,eta)
    358343
    359         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    360         !!!
    361         !!! Purpose: Compute the  constriction of vapor flux by pore ice
    362         !!!
    363         !!! Author: LL
    364         !!!
    365         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    366         implicit none
    367         real,intent(in) :: porefill ! pore filling fraction
    368         real,intent(out) :: eta ! constriction
    369 
    370         !!!
    371 
    372         if (porefill.le.0.) eta = 1.
    373         if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
    374         eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
    375         endif
    376         if (porefill.le.1.) eta = 0.
    377         return
    378         end subroutine
    379 
    380 end module
     344!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     345!!!
     346!!! Purpose: Compute the constriction of vapor flux by pore ice
     347!!!
     348!!! Author: LL
     349!!!
     350!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     351
     352implicit none
     353
     354real, intent(in) :: porefill ! pore filling fraction
     355real, intent(out) :: eta ! constriction
     356
     357if (porefill <= 0.) then
     358    eta = 1.
     359else if (0 < porefill .and. porefill < 1.) then
     360    eta = (1-porefill)**2 ! Hudson et al., JGR, 2009
     361else
     362    eta = 0.
     363endif
     364
     365END SUBROUTINE
     366
     367END MODULE

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

    r3319 r3320  
    3030! Stratum = node of the linked list
    3131type :: stratum
    32     real               :: thickness          ! Layer thickness [m]
    33     real               :: top_elevation      ! Layer top_elevation (top height from the surface) [m]
    34     real               :: co2ice_volfrac    ! CO2 ice volumetric fraction
    35     real               :: h2oice_volfrac    ! H2O ice volumetric fraction
    36     real               :: dust_volfrac       ! Dust volumetric fraction
    37     real               :: air_volfrac        ! Air volumetric fraction inside pores
     32    real                   :: thickness      ! Layer thickness [m]
     33    real                   :: top_elevation  ! Layer top_elevation (top height from the surface) [m]
     34    real                   :: co2ice_volfrac ! CO2 ice volumetric fraction
     35    real                   :: h2oice_volfrac ! H2O ice volumetric fraction
     36    real                   :: dust_volfrac   ! Dust volumetric fraction
     37    real                   :: air_volfrac    ! Air volumetric fraction inside pores
    3838    type(stratum), pointer :: up => null()   ! Upper stratum (next node)
    3939    type(stratum), pointer :: down => null() ! Lower stratum (previous node)
Note: See TracChangeset for help on using the changeset viewer.