Ignore:
Timestamp:
May 2, 2024, 5:06:00 PM (7 months ago)
Author:
jbclement
Message:

PEM:
Correction in the layering algorithm in case of a stratum which is disappearing + few cleanings.
JBC

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

Legend:

Unmodified
Added
Removed

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

    r3320 r3321  
    293293== 30/04/2024 == JBC
    294294Small correction of a bug in the reading of "run_PEM.def" + some cleanings.
     295
     296== 02/05/2024 == JBC
     297Correction in the layering algorithm in case of a stratum which is disappearing + few cleanings.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/launch_pem.sh

    r3319 r3321  
    1818fi
    1919
    20 # Save the current value of LC_NUMERIC and set it to a locale that uses a dot as the decimal separator
    21 OLD_LC_NUMERIC=$LC_NUMERIC
    22 LC_NUMERIC=en_US.UTF-8
    23 
    2420
    2521#####################################################################
    … …  
    5046#####################################################################
    5147
     48
     49# Save the current value of LC_NUMERIC and set it to a locale that uses a dot as the decimal separator
     50OLD_LC_NUMERIC=$LC_NUMERIC
     51LC_NUMERIC=en_US.UTF-8
    5252
    5353#------ Check if files/directories necessary for the script exist ------
    … …  
    167167        mv Xdiurnalave.nc data2reshape_Y${k}.nc
    168168    fi
    169     cp restartfi.nc starts/startfi${iPCM}.nc
     169    cp restartfi.nc starts/restartfi${iPCM}.nc
    170170    mv restartfi.nc startfi.nc
    171171    if [ -f "restart.nc" ]; then
    … …  
    201201    mv diagsoilpem.nc diags/diagsoilpem${iPEM}.nc
    202202fi
    203 cp restartpem.nc starts/startpem${iPEM}.nc
     203cp restartpem.nc starts/restartpem${iPEM}.nc
    204204mv restartpem.nc startpem.nc
    205 cp restartfi.nc starts/startfi_postPEM${iPEM}.nc
     205cp restartfi.nc starts/restartfi_postPEM${iPEM}.nc
    206206mv restartfi.nc startfi.nc
    207207if [ -f "restart.nc" ]; then
    … …  
    245245            mv Xdiurnalave.nc data2reshape_Y${k}.nc
    246246        fi
    247         cp restartfi.nc starts/startfi${iPCM}.nc
     247        cp restartfi.nc starts/restartfi${iPCM}.nc
    248248        mv restartfi.nc startfi.nc
    249249        if [ -f "restart.nc" ]; then
    … …  
    279279        mv diagsoilpem.nc diags/diagsoilpem${iPEM}.nc
    280280    fi
    281     cp restartpem.nc starts/startpem${iPEM}.nc
     281    cp restartpem.nc starts/restartpem${iPEM}.nc
    282282    mv restartpem.nc startpem.nc
    283     cp restartfi.nc starts/startfi_postPEM${iPEM}.nc
     283    cp restartfi.nc starts/restartfi_postPEM${iPEM}.nc
    284284    mv restartfi.nc startfi.nc
    285285    if [ -f "restart.nc" ]; then

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

    r3320 r3321  
    1515real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m]
    1616
    17 !----------------------------------------
     17!-----------------------------------------------------------------------
    1818  contains
    19 !----------------------------------------
     19!-----------------------------------------------------------------------
    2020SUBROUTINE ini_ice_table_porefilling(ngrid,nslope)
    2121
    … …  
    3030END SUBROUTINE ini_ice_table_porefilling
    3131
    32 !----------------------------------------
     32!-----------------------------------------------------------------------
    3333SUBROUTINE end_ice_table_porefilling
    3434
    … …  
    3838if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
    3939
    40  END SUBROUTINE end_ice_table_porefilling
    41 
    42 !----------------------------------------
     40END SUBROUTINE end_ice_table_porefilling
     41
     42!------------------------------------------------------------------------------------------------------
    4343SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
    4444!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    7171!   Locals
    7272! --------
    73 integer ig, islope,isoil,isoilend                              ! loop variables
    74 real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    75 real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
    76 real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
    77 real :: stretch                                                ! stretch factor to improve the convergence of the ice table
    78 real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]
     73integer                       :: ig, islope, isoil, isoilend ! loop variables
     74real, dimension(nsoil_PEM)    :: diff_rho                    ! 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, dimension(ngrid,nslope) :: previous_icetable_depth     ! 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]
    7979!   Code
    8080! ------
    81      previous_icetable_depth(:,:) = ice_table_beg(:,:)
     81previous_icetable_depth = ice_table_beg
    8282do ig = 1,ngrid
    8383    if (watercaptag(ig)) then
    … …  
    9393            if (diff_rho(1) > 0) then ! ice is at the surface
    9494                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
     95                do isoilend = 2,nsoil_PEM - 1
     96                    if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
    9797                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
    9898                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
    … …  
    106106                        ! Now let's find the end of the ice table
    107107                        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
     108                        do isoilend = isoil + 1,nsoil_PEM - 1
    109109                            if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
    110110                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
    … …  
    113113                            endif
    114114                        enddo
    115                     endif !diff_rho(z) <0 & diff_rho(z+1) > 0
     115                    endif ! diff_rho(z) <0 & diff_rho(z+1) > 0
    116116                enddo
    117117            endif ! diff_rho(1) > 0
    … …  
    133133enddo
    134134
    135 END SUBROUTINE
    136 
    137 !----------------------------------------
     135END SUBROUTINE computeice_table_equilibrium
     136
     137!-----------------------------------------------------------------------
    138138SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
    139139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    169169!   Locals
    170170!---------
    171 integer :: ig, islope, ilay, iref    ! loop index
    172 real, dimension(ngrid,nslope) :: rho  ! density of water ice [kg/m^3]
    173 real, dimension(ngrid,nslope) :: Tice ! ice temperature [k]
     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
    175175! ------
    … …  
    187187do ig = 1,ngrid
    188188    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
     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
    190190                          *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    191191    enddo
    192192enddo
    193193
    194 END SUBROUTINE
    195 
    196 !----------------------------------------
    197 
     194END SUBROUTINE compute_massh2o_exchange_ssi
     195
     196!-----------------------------------------------------------------------
    198197SUBROUTINE 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)
    199198!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    245244! 0. Compute constriction over the layer
    246245! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
    247 if (index_IS.lt.0) then
     246if (index_IS < 0) then
    248247    index_tmp = nsoilmx_PEM
    249248    do ilay = 1,nsoilmx_PEM
    … …  
    252251else
    253252    index_tmp = index_IS
    254     do ilay = 1,index_IS-1
     253    do ilay = 1,index_IS - 1
    255254        call constriction(porefill(ilay),eta(ilay))
    256255    enddo
    … …  
    316315endif
    317316if (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)
     317    call  colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,index_geothermal - 1,nsoilmx_PEM,massfillabove)
    319318    index_tmp = -1
    320319    do ilay = index_geothermal,nsoilmx_PEM
    321320        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)
     321        call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
    323322        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
    324323            if (ilay > index_geothermal) then
    … …  
    337336end if
    338337
    339 END SUBROUTINE
    340 
    341 !----------------------------------------
     338END SUBROUTINE find_layering_icetable
     339
     340!-----------------------------------------------------------------------
    342341SUBROUTINE constriction(porefill,eta)
    343342
    … …  
    363362endif
    364363
    365 END SUBROUTINE
    366 
    367 END MODULE
     364END SUBROUTINE constriction
     365
     366END MODULE ice_table_mod

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

    r3320 r3321  
    662662        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
    663663hsubl_dust = r_subl*str%dust_volfrac*str%thickness
    664 ! Update of properties in the sublimating stratum
     664
    665665str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
    666 str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
    667 str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
    668 str%h2oice_volfrac = h_h2oice_old/str%thickness
    669 str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    670 str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
    671 ! Correct the value of top_elevation for the upper strata
    672 current => str%up
    673 do while (associated(current))
    674     current%top_elevation = current%down%top_elevation + current%thickness
    675     current => current%up
    676 enddo
    677 ! Remove the sublimating stratum if there is no ice anymore
    678 if (str%thickness < tol) call rm_stratum(this,str)
     666if (str%thickness < tol) then
     667    ! Remove the sublimating stratum if there is no ice anymore
     668    call rm_stratum(this,str)
     669else
     670    ! Update of properties in the sublimating stratum
     671    str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
     672    str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
     673    str%h2oice_volfrac = h_h2oice_old/str%thickness
     674    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
     675    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     676    ! Correct the value of top_elevation for the upper strata
     677    current => str%up
     678    do while (associated(current))
     679        current%top_elevation = current%down%top_elevation + current%thickness
     680        current => current%up
     681    enddo
     682endif
     683
     684
    679685! New stratum = dust lag
    680686h_lag = hsubl_dust/(1. - dry_lag_porosity)
    … …  
    712718        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
    713719hsubl_dust = r_subl*str%dust_volfrac*str%thickness
    714 ! Update of properties in the sublimating stratum
     720
    715721str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
    716 str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
    717 str%co2ice_volfrac = h_co2ice_old/str%thickness
    718 str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
    719 str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    720 str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
    721 ! Correct the value of top_elevation for the upper strata
    722 current => str%up
    723 do while (associated(current))
    724     current%top_elevation = current%down%top_elevation + current%thickness
    725     current => current%up
    726 enddo
    727 ! Remove the sublimating stratum if there is no ice anymore
    728 if (str%thickness < tol) call rm_stratum(this,str)
     722if (str%thickness < tol) then
     723    ! Remove the sublimating stratum if there is no ice anymore
     724    call rm_stratum(this,str)
     725else
     726    ! Update of properties in the sublimating stratum
     727    str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
     728    str%co2ice_volfrac = h_co2ice_old/str%thickness
     729    str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
     730    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
     731    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     732    ! Correct the value of top_elevation for the upper strata
     733    current => str%up
     734    do while (associated(current))
     735        current%top_elevation = current%down%top_elevation + current%thickness
     736        current => current%up
     737    enddo
     738endif
     739
    729740! New stratum = dust lag
    730741h_lag = hsubl_dust/(1. - dry_lag_porosity)
Note: See TracChangeset for help on using the changeset viewer.