Ignore:
Timestamp:
Mar 7, 2025, 10:01:32 AM (3 months ago)
Author:
jbclement
Message:

PEM:

  • Correction for the update of average pressure at the end.
  • Addition of pore-filled ice information for 'stratum' in the layering algorithm.
  • Deactivating the soil temperature shifting because of bad values given by interpolation.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3667 r3673  
    616616== 03/03/2025 == JBC
    617617Adjusting the anticipation time for more security regarding the job time limit detection.
     618
     619== 07/03/2025 == JBC
     620- Correction for the update of average pressure at the end.
     621- Addition of pore-filled ice information for 'stratum' in the layering algorithm.
     622- Deactivating the soil temperature shifting because of bad values given by interpolation.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

    r3666 r3673  
    6767    stratif_h2oice = np.zeros((ngrid,nfile,nslope,max_nb_str))
    6868    stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str))
    69     stratif_air = np.zeros((ngrid,nfile,nslope,max_nb_str))
     69    stratif_pore = np.zeros((ngrid,nfile,nslope,max_nb_str))
    7070    for var_name in datasets[0].variables:
    7171        if 'top_elevation' in var_name:
    … …  
    9797                            stratif_dust[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    9898            print(f"Processed variable: {var_name}")
    99         elif 'air_volfrac' in var_name:
     99        elif 'pore_volfrac' in var_name:
    100100            for i in range(0,ngrid):
    101101                for j in range(0,nfile):
    102102                    for k in range(0,nslope):
    103103                        if f'slope{k + 1:02d}' in var_name:
    104                             stratif_air[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
     104                            stratif_pore[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    105105            print(f"Processed variable: {var_name}")
    106106
    … …  
    109109        ds.close()
    110110
    111     stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_air]
     111    stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_pore]
    112112    return stratif_data, max_top_elevation, longitude, latitude
    113113
    … …  
    158158
    159159### Figures plotting
    160 subtitle = ['CO2 ice','H2O ice','Dust','Air']
     160subtitle = ['CO2 ice','H2O ice','Dust','Pore']
    161161cmap = plt.get_cmap('viridis').copy()
    162162cmap.set_under('white')

  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py

    r3458 r3673  
    4343stratif_h2oice_volfrac = []
    4444stratif_dust_volfrac = []
    45 stratif_air_volfrac = []
     45stratif_pore_volfrac = []
    4646for i in range(1,nslope + 1):
    4747    stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:])
    … …  
    5050    stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:])
    5151    stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:])
    52     stratif_air_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_air_volfrac'][:])
     52    stratif_pore_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_pore_volfrac'][:])
    5353
    5454### Display the data
    5555igrid = igrid - 1
    5656islope = islope - 1
    57 labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air'
     57labels = 'CO2 ice', 'H2O ice', 'Dust', 'Pore'
    5858contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1])
    5959height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1)
    … …  
    6464    contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]
    6565    contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]
    66     contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid]
     66    contents[3,1 + i] = stratif_pore_volfrac[islope][0,i,igrid]
    6767contents[:,0] = contents[:,1]
    6868

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

    r3666 r3673  
    3333! Stratum = node of the linked list
    3434type :: stratum
    35     real                   :: thickness      ! Layer thickness [m]
    36     real                   :: top_elevation  ! Layer top_elevation (top height from the surface) [m]
    37     real                   :: co2ice_volfrac ! CO2 ice volumetric fraction
    38     real                   :: h2oice_volfrac ! H2O ice volumetric fraction
    39     real                   :: dust_volfrac   ! Dust volumetric fraction
    40     real                   :: air_volfrac    ! Air volumetric fraction inside pores
    41     type(stratum), pointer :: up => null()   ! Upper stratum (next node)
    42     type(stratum), pointer :: down => null() ! Lower stratum (previous node)
     35    real                   :: thickness       ! Layer thickness [m]
     36    real                   :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
     37    real                   :: co2ice_volfrac  ! CO2 ice volumetric fraction
     38    real                   :: h2oice_volfrac  ! H2O ice volumetric fraction
     39    real                   :: dust_volfrac    ! Dust volumetric fraction
     40    real                   :: pore_volfrac    ! Pore volumetric fraction = pores
     41    real                   :: poreice_volfrac ! Pore-filled ice volumetric fraction
     42    type(stratum), pointer :: up => null()    ! Upper stratum (next node)
     43    type(stratum), pointer :: down => null()  ! Lower stratum (previous node)
    4344end type stratum
    4445
    … …  
    8586
    8687!---- Code
    87 write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
    88 write(*,'(a20,es13.6)') 'Top elevation  : ', this%top_elevation
    89 write(*,'(a20,es13.6)') 'CO2 ice (m3/m3): ', this%co2ice_volfrac
    90 write(*,'(a20,es13.6)') 'H2O ice (m3/m3): ', this%h2oice_volfrac
    91 write(*,'(a20,es13.6)') 'Dust (m3/m3)   : ', this%dust_volfrac
    92 write(*,'(a20,es13.6)') 'Air (m3/m3)    : ', this%air_volfrac
     88write(*,'(a,es13.6)') 'Thickness              : ', this%thickness
     89write(*,'(a,es13.6)') 'Top elevation          : ', this%top_elevation
     90write(*,'(a,es13.6)') 'CO2 ice (m3/m3)        : ', this%co2ice_volfrac
     91write(*,'(a,es13.6)') 'H2O ice (m3/m3)        : ', this%h2oice_volfrac
     92write(*,'(a,es13.6)') 'Dust (m3/m3)           : ', this%dust_volfrac
     93write(*,'(a,es13.6)') 'Porosity (m3/m3)       : ', this%pore_volfrac
     94write(*,'(a,es13.6)') 'Pore-filled ice (m3/m3): ', this%poreice_volfrac
    9395
    9496END SUBROUTINE print_stratum
    … …  
    9698!=======================================================================
    9799! To add a stratum to the top of a layering
    98 SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
     100SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
    99101
    100102implicit none
    … …  
    102104!---- Arguments
    103105type(layering), intent(inout)  :: this
    104 real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
     106real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
    105107
    106108!---- Local variables
    … …  
    116118str%h2oice_volfrac = 0.
    117119str%dust_volfrac = 0.
    118 str%air_volfrac = 1.
     120str%pore_volfrac = 1.
     121str%poreice_volfrac = 0.
    119122if (present(thickness)) str%thickness = thickness
    120123if (present(top_elevation)) str%top_elevation = top_elevation
    … …  
    122125if (present(h2oice)) str%h2oice_volfrac = h2oice
    123126if (present(dust)) str%dust_volfrac = dust
    124 if (present(air)) str%air_volfrac = air
     127if (present(pore)) str%pore_volfrac = pore
     128if (present(poreice)) str%poreice_volfrac = poreice
    125129
    126130! Verification of volume fraction
    127 if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) then
     131if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
    128132    call print_stratum(str)
    129133    error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
    … …  
    147151!=======================================================================
    148152! To insert a stratum after another one in a layering
    149 SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,air)
     153SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
    150154
    151155implicit none
    … …  
    154158type(layering),         intent(inout) :: this
    155159type(stratum), pointer, intent(inout) :: str_prev
    156 real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
     160real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
    157161
    158162!---- Local variables
    … …  
    161165!---- Code
    162166if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
    163     call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
     167    call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
    164168else
    165169    ! Creation of the stratum
    … …  
    171175    str%h2oice_volfrac = 0.
    172176    str%dust_volfrac = 0.
    173     str%air_volfrac = 1.
     177    str%pore_volfrac = 1.
     178    str%poreice_volfrac = 0.
    174179    if (present(thickness)) str%thickness = thickness
    175180    if (present(top_elevation)) str%top_elevation = top_elevation
    … …  
    177182    if (present(h2oice)) str%h2oice_volfrac = h2oice
    178183    if (present(dust)) str%dust_volfrac = dust
    179     if (present(air)) str%air_volfrac = air
     184    if (present(pore)) str%pore_volfrac = pore
     185    if (present(poreice)) str%poreice_volfrac = poreice
    180186
    181187    ! Verification of volume fraction
    182     if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) then
     188    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
    183189        call print_stratum(str)
    184190        error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
    … …  
    326332i = this%nb_str
    327333do while (associated(current))
    328     write(*,'(a8,i4)') 'Stratum ', i
     334    write(*,'(a,i4)') 'Stratum ', i
    329335    call print_stratum(current)
    330336    write(*,*)
    … …  
    386392            stratif_array(ig,islope,k,4) = current%h2oice_volfrac
    387393            stratif_array(ig,islope,k,5) = current%dust_volfrac
    388             stratif_array(ig,islope,k,6) = current%air_volfrac
     394            stratif_array(ig,islope,k,6) = current%pore_volfrac
     395            stratif_array(ig,islope,k,7) = current%poreice_volfrac
    389396            current => current%up
    390397            k = k + 1
    … …  
    412419do islope = 1,nslope
    413420    do ig = 1,ngrid
    414         stratif(ig,islope)%bottom%thickness      = stratif_array(ig,islope,1,1)
    415         stratif(ig,islope)%bottom%top_elevation  = stratif_array(ig,islope,1,2)
    416         stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3)
    417         stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4)
    418         stratif(ig,islope)%bottom%dust_volfrac   = stratif_array(ig,islope,1,5)
    419         stratif(ig,islope)%bottom%air_volfrac    = stratif_array(ig,islope,1,6)
     421        stratif(ig,islope)%bottom%thickness       = stratif_array(ig,islope,1,1)
     422        stratif(ig,islope)%bottom%top_elevation   = stratif_array(ig,islope,1,2)
     423        stratif(ig,islope)%bottom%co2ice_volfrac  = stratif_array(ig,islope,1,3)
     424        stratif(ig,islope)%bottom%h2oice_volfrac  = stratif_array(ig,islope,1,4)
     425        stratif(ig,islope)%bottom%dust_volfrac    = stratif_array(ig,islope,1,5)
     426        stratif(ig,islope)%bottom%pore_volfrac    = stratif_array(ig,islope,1,6)
     427        stratif(ig,islope)%bottom%poreice_volfrac = stratif_array(ig,islope,1,7)
    420428        do k = 2,size(stratif_array,3)
    421429            if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
    422             call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
     430            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6),stratif_array(ig,islope,k,7))
    423431        enddo
    424432    enddo
    … …  
    490498    do while (h2lift_tot > 0. .and. associated(current1))
    491499        ! Is the considered stratum a dust lag?
    492         if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we can lift dust
     500        if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we can lift dust
    493501            h_dust_old = current1%dust_volfrac*current1%thickness
    494502            ! How much dust can be lifted in the considered dust lag?
    … …  
    524532        if (thickness > 0.) then ! Only if the stratum is building up
    525533             if (new_str) then
    526                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
     534                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity,0.)
    527535                 new_str = .false.
    528536             else
    … …  
    540548        if (thickness > 0.) then ! Only if the stratum is building up
    541549             if (new_str) then
    542                  call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
     550                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness),0.)
    543551                 new_str = .false.
    544552             else
    … …  
    563571            if (associated(current1%up)) then ! If there is an upper stratum
    564572                currentloc => current1%up
    565                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
     573                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    566574                    h_toplag = h_toplag + currentloc%thickness
    567575                    if (.not. associated(currentloc%up)) exit
    … …  
    578586            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    579587                ! Is the considered stratum a dust lag?
    580                 if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
     588                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    581589                    current1 => current1%down
    582590                    new_lag1 = .true.
    … …  
    600608        if (thickness > 0.) then ! Only if the stratum is building up
    601609             if (new_str) then
    602                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
     610                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness),0.)
    603611                 new_str = .false.
    604612             else
    … …  
    623631            if (associated(current1%up)) then ! If there is an upper stratum
    624632                currentloc => current1%up
    625                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
     633                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    626634                    h_toplag = h_toplag + currentloc%thickness
    627635                    if (.not. associated(currentloc%up)) exit
    … …  
    638646            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    639647                ! Is the considered stratum a dust lag?
    640                 if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
     648                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    641649                    current1 => current1%down
    642650                    new_lag1 = .true.
    … …  
    667675            if (associated(current2%up)) then ! If there is an upper stratum
    668676                currentloc => current2%up
    669                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
     677                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    670678                    h_toplag = h_toplag + currentloc%thickness
    671679                    if (.not. associated(currentloc%up)) exit
    … …  
    682690            else if (h_h2oice_old < eps) then ! There is nothing to sublimate
    683691                ! Is the considered stratum a dust lag?
    684                 if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
     692                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    685693                    current1 => current1%down
    686694                    new_lag1 = .true.
    … …  
    743751    str%h2oice_volfrac = h_h2oice_old/str%thickness
    744752    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    745     str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     753    str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     754    str%poreice_volfrac = 0.
    746755    ! Correct the value of top_elevation for the upper strata
    747756    current => str%up
    … …  
    801810    str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
    802811    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    803     str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     812    str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
     813    str%poreice_volfrac = 0.
    804814    ! Correct the value of top_elevation for the upper strata
    805815    current => str%up
    … …  
    838848!---- Code
    839849! Update of properties in the eroding dust lag
    840 str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
    841 str%top_elevation = str%top_elevation - h2lift/(1. - str%air_volfrac)
     850str%thickness = str%thickness - h2lift/(1. - str%pore_volfrac)
     851str%top_elevation = str%top_elevation - h2lift/(1. - str%pore_volfrac)
    842852! Remove the eroding dust lag if there is no dust anymore
    843853if (str%thickness < tol) call rm_stratum(this,str)

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

    r3667 r3673  
    145145real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
    146146real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
    147 real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Averaged surface pressure
     147real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Average surface pressure
    148148real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
    149149real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
    … …  
    206206
    207207! Variables for surface and soil
    208 real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Averaged surface temperature [K]
     208real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Average surface temperature [K]
    209209real, dimension(:,:),     allocatable :: tsurf_dev                        ! Grid points x Slope field: Surface temperature deviation [K]
    210 real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Averaged surface temperature of first call of the PCM [K]
    211 real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Averaged Soil Temperature [K]
    212 real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Averaged Soil Temperature at the previous time step [K]
     210real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
     211real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Average Soil Temperature [K]
     212real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
    213213real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
    214214real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
    … …  
    216216real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old         ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
    217217real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
    218 real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Averaged water surface density [kg/m^3]
     218real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Average water surface density [kg/m^3]
    219219real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
    220 real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
     220real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
    221221real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    222222real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    … …  
    374374call hostnm(hostname) ! Machine/station name
    375375write(*,*)
    376 write(*,*) '######### PEM information #########'
     376write(*,*) '********* PEM information *********'
    377377write(*,*) '+ User     : '//trim(logname)
    378378write(*,*) '+ Machine  : '//trim(hostname)
    … …  
    406406!------------------------
    407407write(*,*)
    408 write(*,*) '######### PEM initialization #########'
     408write(*,*) '********* PEM initialization *********'
    409409write(*,*) '> Reading "run.def" (PEM)'
    410410#ifndef CPP_1D
    … …  
    638638ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
    639639ps_avg_global_new = ps_avg_global_ini
    640 write(*,*) "Total surface of the planet      =", total_surface
    641 write(*,*) "Initial global averaged pressure =", ps_avg_global_ini
     640write(*,*) "Total surface of the planet     =", total_surface
     641write(*,*) "Initial global average pressure =", ps_avg_global_ini
    642642
    643643!------------------------
    … …  
    732732!------------------------
    733733write(*,*)
    734 write(*,*) '######### PEM cycle #########'
     734write(*,*) '********* PEM cycle *********'
    735735i_myear_leg = 0
    736736stopPEM = 0
    … …  
    750750do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
    751751! II.a.1. Compute updated global pressure
    752     write(*,'(a,f10.2)') ' #### Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
     752    write(*,'(a,f10.2)') ' **** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
    753753    write(*,*) "> Updating the surface pressure"
    754754    ps_avg_global_old = ps_avg_global_new
    … …  
    775775        enddo
    776776    enddo
    777     do ig = 1,ngrid
    778         ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old
    779     enddo
     777    ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
    780778    write(*,*) "> Updating the pressure levels timeseries for the new pressure"
    781779    allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
    … …  
    893891    if (soil_pem) then
    894892! II_d.2 Shifting soil temperature to surface
    895         call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
     893        !call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
    896894
    897895! II_d.3 Update soil temperature
    … …  
    939937            enddo
    940938            deallocate(porefill)
    941             call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
     939            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    942940        endif
    943941        deallocate(icetable_thickness_old,ice_porefilling_old)
    … …  
    10401038        select case (stopPEM)
    10411039            case(1)
    1042                 write(*,'(a,i0,a)') " #### STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
     1040                write(*,'(a,i0,a)') " **** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
    10431041            case(2)
    1044                 write(*,'(a,i0,a)') " #### STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
     1042                write(*,'(a,i0,a)') " **** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
    10451043            case(3)
    1046                 write(*,'(a,i0,a)') " #### STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
     1044                write(*,'(a,i0,a)') " **** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
    10471045            case(4)
    1048                 write(*,'(a,i0,a)') " #### STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
     1046                write(*,'(a,i0,a)') " **** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
    10491047            case(5)
    1050                 write(*,'(a,i0)') " #### STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
     1048                write(*,'(a,i0)') " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
    10511049            case(6)
    1052                 write(*,'(a,i0)') " #### STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
     1050                write(*,'(a,i0)') " **** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
    10531051            case(7)
    1054                 write(*,'(a,i0)') " #### STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
     1052                write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
    10551053            case(8)
    1056                 write(*,'(a,i0)') " #### STOPPING because the layering algorithm met an hasty end: ", stopPEM
     1054                write(*,'(a,i0)') " **** STOPPING because the layering algorithm met an hasty end: ", stopPEM
    10571055            case default
    1058                 write(*,'(a,i0)') " #### STOPPING with unknown stopping criterion code: ", stopPEM
     1056                write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM
    10591057        end select
    10601058        exit
    10611059    else
    1062         write(*,'(a,f10.2,a)') ' #### The chained simulation has run for ',i_myear,' Martian years.'
    1063         write(*,*) '#### The PEM can continue!'
    1064         write(*,*) '####'
     1060        write(*,'(a,f10.2,a)') ' **** The chained simulation has run for ',i_myear,' Martian years.'
     1061        write(*,*) '**** The PEM can continue!'
     1062        write(*,*) '****'
    10651063    endif
    10661064enddo ! big time iteration loop of the pem
    … …  
    10811079! III_a.1 Ice update for start file
    10821080write(*,*)
    1083 write(*,*) '######### PEM finalization #########'
     1081write(*,*) '********* PEM finalization *********'
    10841082write(*,*) '> Reconstructing perennial ice and frost for the PCM'
    10851083watercap = 0.
    … …  
    11331131write(*,*) '> Reconstructing the pressure for the PCM'
    11341132allocate(ps_start(ngrid))
    1135 ! The pressure deviation is rescaled to avoid disproportionate oscillations in case of huge average pressure drop
    1136 ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
     1133ps_start = (ps_avg + ps_dev)*ps_avg_global_new/ps_avg_global_ini ! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
    11371134deallocate(ps_avg,ps_dev)
    11381135
    … …  
    12881285
    12891286write(*,*)
    1290 write(*,*) '######### PEM finalization #########'
     1287write(*,*) '********* PEM finalization *********'
    12911288write(*,'(a,f10.2,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
    12921289write(*,'(a,f10.2,a,f10.2,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years."
    12931290write(*,'(a,f10.2,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years."
    12941291write(*,*) "+ PEM: so far, so good!"
    1295 write(*,*) '####################################'
     1292write(*,*) '************************************'
    12961293
    12971294if (layering_algo) then

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

    r3666 r3673  
    6565real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
    6666! local
    67 real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem           ! soil temperature saved in the start [K]
    68 real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM              ! soil thermal inertia saved in the start [SI]
    69 logical                                 :: found, found2            ! check if variables are found in the start
    70 integer                                 :: iloop, ig, islope, isoil ! index for loops
    71 real                                    :: delta                    ! Depth of the interface regolith-breccia, breccia -bedrock [m]
    72 character(2)                            :: num                      ! intermediate string to read PEM start sloped variables
    73 logical                                 :: startpem_file            ! boolean to check if we read the startfile or not
    74 real, dimension(:,:,:,:), allocatable   :: stratif_array            ! Array for stratification (layerings)
     67real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem    ! soil temperature saved in the start [K]
     68real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM       ! soil thermal inertia saved in the start [SI]
     69logical                                 :: found, found2     ! check if variables are found in the start
     70integer                                 :: iloop, ig, islope ! index for loops
     71real                                    :: delta             ! Depth of the interface regolith-breccia, breccia -bedrock [m]
     72character(2)                            :: num               ! intermediate string to read PEM start sloped variables
     73logical                                 :: startpem_file     ! boolean to check if we read the startfile or not
     74real, dimension(:,:,:,:), allocatable   :: stratif_array     ! Array for stratification (layerings)
    7575
    7676#ifdef CPP_STD
    … …  
    143143            nb_str_max = int(inquire_dimension_length('nb_str_max'))
    144144        endif
    145         allocate(stratif_array(ngrid,nslope,nb_str_max,6))
     145        allocate(stratif_array(ngrid,nslope,nb_str_max,7))
    146146        stratif_array = 0.
    147147        do islope = 1,nslope
    … …  
    156156            call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
    157157            found = found .or. found2
    158             call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
     158            call get_field('stratif_slope'//num//'_pore_volfrac',stratif_array(:,islope,:,6),found2)
     159            found = found .or. found2
     160            call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,7),found2)
    159161            found = found .or. found2
    160162            if (.not. found) then
    … …  
    267269                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
    268270
    269                 do iloop = nsoil_PCM + 1,nsoil_PEM
    270                     tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
    271                 enddo
     271                tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) = tsoil_startpem(:,nsoil_PCM + 1:nsoil_PEM,islope)
    272272            endif !found
    273273
    274             do isoil = nsoil_PCM + 1,nsoil_PEM
    275                 watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    276             enddo
     274            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    277275        enddo ! islope
    278276        write(*,*) 'PEMETAT0: TSOIL done'
    … …  
    431429
    432430! First raw initialization
    433             do isoil = nsoil_PCM + 1,nsoil_PEM
    434                 do ig = 1,ngrid
    435                     watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    436                 enddo
    437             enddo
     431            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    438432        enddo !islope
    439433        write(*,*) 'PEMETAT0: TSOIL done'

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

    r3591 r3673  
    9999
    100100if (layering_algo) then
    101     allocate(stratif_array(ngrid,nslope,nb_str_max,6))
     101    allocate(stratif_array(ngrid,nslope,nb_str_max,7))
    102102    call stratif2array(stratif,ngrid,nslope,stratif_array)
    103103    do islope = 1,nslope
    … …  
    108108        call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
    109109        call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
    110         call put_field('stratif_slope'//num//'_air_volfrac','Layering air volume fraction',stratif_array(:,islope,:,6),Year)
     110        call put_field('stratif_slope'//num//'_pore_volfrac','Layering pore volume fraction',stratif_array(:,islope,:,6),Year)
     111        call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,7),Year)
    111112    enddo
    112113    deallocate(stratif_array)
Note: See TracChangeset for help on using the changeset viewer.