[3076] | 1 | MODULE pemetat0_mod |
---|
[3019] | 2 | |
---|
[3076] | 3 | implicit none |
---|
| 4 | |
---|
| 5 | !======================================================================= |
---|
| 6 | contains |
---|
| 7 | !======================================================================= |
---|
| 8 | |
---|
[3525] | 9 | SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, & |
---|
| 10 | ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, & |
---|
| 11 | global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & |
---|
[3297] | 12 | m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif) |
---|
[3076] | 13 | |
---|
[3297] | 14 | use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length |
---|
[3161] | 15 | use comsoil_h_PEM, only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock |
---|
[3149] | 16 | use comsoil_h, only: volcapa, inertiedat |
---|
[3076] | 17 | use adsorption_mod, only: regolith_adsorption, adsorption_pem |
---|
[3525] | 18 | use ice_table_mod, only: computeice_table_equilibrium, icetable_depth, icetable_thickness, icetable_equilibrium, icetable_dynamic |
---|
[3076] | 19 | use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock |
---|
| 20 | use soil_thermalproperties_mod, only: update_soil_thermalproperties |
---|
| 21 | use tracer_mod, only: mmol, igcm_h2o_vap ! tracer names and molar masses |
---|
| 22 | use abort_pem_mod, only: abort_pem |
---|
[3178] | 23 | use compute_soiltemp_mod, only: ini_tsoil_pem, compute_tsoil_pem |
---|
[3161] | 24 | use comslope_mod, only: def_slope_mean, subslope_dist |
---|
[3319] | 25 | use layering_mod, only: layering, array2stratif, nb_str_max, layering_algo |
---|
[3076] | 26 | |
---|
[3019] | 27 | #ifndef CPP_STD |
---|
[3161] | 28 | use comcstfi_h, only: r, mugaz, pi |
---|
| 29 | use surfdat_h, only: watercaptag, watercap, perennial_co2ice |
---|
[3019] | 30 | #else |
---|
[3161] | 31 | use comcstfi_mod, only: r, mugaz, pi |
---|
[3019] | 32 | #endif |
---|
| 33 | |
---|
[3076] | 34 | implicit none |
---|
[3019] | 35 | |
---|
[3076] | 36 | include "callkeys.h" |
---|
[2835] | 37 | |
---|
[3297] | 38 | character(*), intent(in) :: filename ! name of the startfi_PEM.nc |
---|
[3076] | 39 | integer, intent(in) :: ngrid ! # of physical grid points |
---|
[3149] | 40 | integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM |
---|
[3076] | 41 | integer, intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM |
---|
| 42 | integer, intent(in) :: nslope ! # of sub-grid slopes |
---|
| 43 | integer, intent(in) :: timelen ! # time samples |
---|
| 44 | real, intent(in) :: timestep ! time step [s] |
---|
[3149] | 45 | real, intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa] |
---|
| 46 | real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K] |
---|
[3532] | 47 | real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K] |
---|
[3076] | 48 | real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg] |
---|
| 49 | real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg] |
---|
| 50 | real, dimension(ngrid,timelen), intent(in) :: ps_inst ! surface pressure [Pa] |
---|
[3498] | 51 | real, dimension(ngrid,nslope), intent(in) :: d_h2oice ! tendencies on h2o ice |
---|
| 52 | real, dimension(ngrid,nslope), intent(in) :: d_co2ice ! tendencies on co2 ice |
---|
[3149] | 53 | real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3) |
---|
[2794] | 54 | ! outputs |
---|
[3076] | 55 | real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2] |
---|
| 56 | real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2] |
---|
[3149] | 57 | real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! h2o ice amount [kg/m^2] |
---|
| 58 | real, dimension(ngrid,nslope), intent(out) :: co2_ice ! co2 ice amount [kg/m^2] |
---|
[3297] | 59 | type(layering), dimension(ngrid,nslope), intent(inout) :: stratif ! stratification (layerings) |
---|
[3076] | 60 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI] |
---|
| 61 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K] |
---|
[3493] | 62 | real, dimension(ngrid,nslope), intent(inout) :: ice_table_depth ! Ice table depth [m] |
---|
[3076] | 63 | real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m] |
---|
[3493] | 64 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3] |
---|
[3076] | 65 | real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K] |
---|
| 66 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2] |
---|
| 67 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2] |
---|
[3149] | 68 | real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! surface water ice density, yearly averaged (kg/m^3) |
---|
[2794] | 69 | ! local |
---|
[3297] | 70 | real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM ! soil temperature saved in the start [K] |
---|
| 71 | real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! soil thermal inertia saved in the start [SI] |
---|
| 72 | logical :: found, found2 ! check if variables are found in the start |
---|
| 73 | integer :: iloop, ig, islope, it, isoil, k ! index for loops |
---|
| 74 | real :: kcond ! Thermal conductivity, intermediate variable [SI] |
---|
| 75 | real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m] |
---|
| 76 | character(2) :: num ! intermediate string to read PEM start sloped variables |
---|
| 77 | real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1 ! intermediate soil temperature during yr 1 [K] |
---|
| 78 | real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2 ! intermediate soil temperature during yr 2 [K] |
---|
| 79 | logical :: startpem_file ! boolean to check if we read the startfile or not |
---|
| 80 | real, dimension(:,:,:,:), allocatable :: stratif_array ! Array for stratification (layerings) |
---|
[2835] | 81 | |
---|
[3076] | 82 | #ifdef CPP_STD |
---|
[3149] | 83 | logical, dimension(ngrid) :: watercaptag |
---|
[3159] | 84 | watercaptag = .false. |
---|
[2985] | 85 | #endif |
---|
| 86 | |
---|
[2794] | 87 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 88 | !!! |
---|
| 89 | !!! Purpose: read start_pem. Need a specific iostart_PEM |
---|
| 90 | !!! |
---|
[2888] | 91 | !!! Order: 0. Previous year of the PEM run |
---|
[3149] | 92 | !!! Ice initialization |
---|
[3076] | 93 | !!! 1. Thermal Inertia |
---|
[2888] | 94 | !!! 2. Soil Temperature |
---|
| 95 | !!! 3. Ice table |
---|
| 96 | !!! 4. Mass of CO2 & H2O adsorbed |
---|
[2794] | 97 | !!! |
---|
| 98 | !!! /!\ This order must be respected ! |
---|
| 99 | !!! Author: LL |
---|
| 100 | !!! |
---|
| 101 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2779] | 102 | |
---|
[3170] | 103 | !0.1 Check if the start_PEM exist. |
---|
[2863] | 104 | |
---|
[3143] | 105 | inquire(file = filename,exist = startpem_file) |
---|
[2863] | 106 | |
---|
[3498] | 107 | write(*,*)'Is there a "startpem" file?',startpem_file |
---|
[2863] | 108 | |
---|
[3170] | 109 | !0.2 Set to default values |
---|
[3493] | 110 | ice_table_depth = -1. ! by default, no ice table |
---|
[3170] | 111 | ice_table_thickness = -1. |
---|
[2794] | 112 | !1. Run |
---|
| 113 | if (startpem_file) then |
---|
[3214] | 114 | write(*,*)'There is a "startpem.nc"!' |
---|
[3076] | 115 | ! open pem initial state file: |
---|
| 116 | call open_startphy(filename) |
---|
[2835] | 117 | |
---|
[3149] | 118 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[3430] | 119 | #ifndef CPP_STD |
---|
| 120 | ! h2o ice |
---|
| 121 | h2o_ice = 0. |
---|
| 122 | call get_field("h2o_ice",h2o_ice,found) |
---|
| 123 | if (.not. found) then |
---|
| 124 | write(*,*)'Pemetat0: failed loading <h2o_ice>' |
---|
| 125 | write(*,*)'will reconstruct the values from watercaptag' |
---|
| 126 | write(*,*)'with default value ''ini_huge_h2oice''' |
---|
| 127 | do ig = 1,ngrid |
---|
| 128 | if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice |
---|
| 129 | enddo |
---|
| 130 | else |
---|
| 131 | ! The variations of infinite reservoirs during the PCM years are taken into account |
---|
| 132 | h2o_ice = h2o_ice + watercap |
---|
| 133 | endif |
---|
| 134 | |
---|
| 135 | ! co2 ice |
---|
| 136 | co2_ice = perennial_co2ice |
---|
| 137 | #endif |
---|
| 138 | |
---|
| 139 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[3297] | 140 | ! Stratification (layerings) |
---|
[3430] | 141 | nb_str_max = 1 |
---|
[3319] | 142 | if (layering_algo) then |
---|
| 143 | found = inquire_dimension("nb_str_max") |
---|
| 144 | if (.not. found) then |
---|
| 145 | write(*,*) 'Pemetat0: failed loading <nb_str_max>!' |
---|
| 146 | write(*,*) '''nb_str_max'' is set to 3!' |
---|
| 147 | else |
---|
| 148 | nb_str_max = int(inquire_dimension_length('nb_str_max')) |
---|
| 149 | endif |
---|
| 150 | allocate(stratif_array(ngrid,nslope,nb_str_max,6)) |
---|
| 151 | stratif_array = 0. |
---|
| 152 | do islope = 1,nslope |
---|
| 153 | write(num,'(i2.2)') islope |
---|
| 154 | call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found) |
---|
| 155 | call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2) |
---|
| 156 | found = found .or. found2 |
---|
| 157 | call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2) |
---|
| 158 | found = found .or. found2 |
---|
| 159 | call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2) |
---|
| 160 | found = found .or. found2 |
---|
| 161 | call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2) |
---|
| 162 | found = found .or. found2 |
---|
| 163 | call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2) |
---|
| 164 | found = found .or. found2 |
---|
| 165 | if (.not. found) then |
---|
| 166 | write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>' |
---|
| 167 | endif ! not found |
---|
| 168 | enddo ! islope |
---|
| 169 | call array2stratif(stratif_array,ngrid,nslope,stratif) |
---|
| 170 | deallocate(stratif_array) |
---|
[3297] | 171 | endif |
---|
| 172 | |
---|
[3039] | 173 | if (soil_pem) then |
---|
[3076] | 174 | |
---|
[3149] | 175 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2794] | 176 | !1. Thermal Inertia |
---|
| 177 | ! a. General case |
---|
[3149] | 178 | do islope = 1,nslope |
---|
| 179 | write(num,'(i2.2)') islope |
---|
| 180 | call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found) |
---|
| 181 | if (.not. found) then |
---|
| 182 | write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>' |
---|
| 183 | write(*,*)'will reconstruct the values of TI_PEM' |
---|
[2779] | 184 | |
---|
[3149] | 185 | do ig = 1,ngrid |
---|
| 186 | if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then |
---|
| 187 | !!! transition |
---|
| 188 | delta = depth_breccia |
---|
[3430] | 189 | TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & |
---|
| 190 | (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + & |
---|
| 191 | ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) |
---|
| 192 | do iloop=index_breccia + 2,index_bedrock |
---|
[3149] | 193 | TI_PEM(ig,iloop,islope) = TI_breccia |
---|
| 194 | enddo |
---|
[3159] | 195 | else ! we keep the high TI values |
---|
[3149] | 196 | do iloop = index_breccia + 1,index_bedrock |
---|
| 197 | TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope) |
---|
| 198 | enddo |
---|
| 199 | endif ! TI PEM and breccia comparison |
---|
| 200 | !! transition |
---|
| 201 | delta = depth_bedrock |
---|
| 202 | TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & |
---|
| 203 | (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + & |
---|
| 204 | ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2)))) |
---|
| 205 | do iloop = index_bedrock + 2,nsoil_PEM |
---|
| 206 | TI_PEM(ig,iloop,islope) = TI_bedrock |
---|
[3076] | 207 | enddo |
---|
| 208 | enddo |
---|
[3149] | 209 | else ! found |
---|
| 210 | do iloop = nsoil_PCM + 1,nsoil_PEM |
---|
| 211 | TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope) ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here. |
---|
| 212 | enddo |
---|
| 213 | endif ! not found |
---|
| 214 | enddo ! islope |
---|
[2779] | 215 | |
---|
[3149] | 216 | write(*,*) 'PEMETAT0: THERMAL INERTIA done' |
---|
[2779] | 217 | |
---|
[2794] | 218 | ! b. Special case for inertiedat, inertiedat_PEM |
---|
[3149] | 219 | call get_field("inertiedat_PEM",inertiedat_PEM,found) |
---|
[3456] | 220 | if (.not. found) then |
---|
[3149] | 221 | do iloop = 1,nsoil_PCM |
---|
| 222 | inertiedat_PEM(:,iloop) = inertiedat(:,iloop) |
---|
| 223 | enddo |
---|
| 224 | !!! zone de transition |
---|
| 225 | delta = depth_breccia |
---|
| 226 | do ig = 1,ngrid |
---|
[3430] | 227 | if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then |
---|
[3493] | 228 | inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & |
---|
| 229 | (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & |
---|
[3149] | 230 | ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) |
---|
[2779] | 231 | |
---|
[3430] | 232 | do iloop = index_breccia + 2,index_bedrock |
---|
[3149] | 233 | inertiedat_PEM(ig,iloop) = TI_breccia |
---|
| 234 | enddo |
---|
| 235 | else |
---|
[3430] | 236 | do iloop=index_breccia + 1,index_bedrock |
---|
[3149] | 237 | inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM) |
---|
| 238 | enddo |
---|
| 239 | endif ! comparison ti breccia |
---|
| 240 | enddo ! ig |
---|
[2779] | 241 | |
---|
[3149] | 242 | !!! zone de transition |
---|
| 243 | delta = depth_bedrock |
---|
| 244 | do ig = 1,ngrid |
---|
[3430] | 245 | inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ & |
---|
| 246 | (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ & |
---|
| 247 | ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2)))) |
---|
[3149] | 248 | enddo |
---|
[2779] | 249 | |
---|
[3149] | 250 | do iloop = index_bedrock + 2, nsoil_PEM |
---|
| 251 | do ig = 1,ngrid |
---|
| 252 | inertiedat_PEM(ig,iloop) = TI_bedrock |
---|
| 253 | enddo |
---|
| 254 | enddo |
---|
| 255 | endif ! not found |
---|
[2779] | 256 | |
---|
[3076] | 257 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2794] | 258 | !2. Soil Temperature |
---|
[3149] | 259 | do islope=1,nslope |
---|
| 260 | write(num,fmt='(i2.2)') islope |
---|
| 261 | call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found) |
---|
[3456] | 262 | if (.not. found) then |
---|
[3149] | 263 | write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>' |
---|
| 264 | write(*,*)'will reconstruct the values of Tsoil' |
---|
| 265 | ! do ig = 1,ngrid |
---|
| 266 | ! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa |
---|
| 267 | ! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1)) |
---|
| 268 | ! do iloop=index_breccia+2,index_bedrock |
---|
| 269 | ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa |
---|
| 270 | ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia)) |
---|
| 271 | ! enddo |
---|
| 272 | ! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa |
---|
| 273 | ! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) |
---|
[2895] | 274 | ! |
---|
[3149] | 275 | ! do iloop=index_bedrock+2,nsoil_PEM |
---|
| 276 | ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa |
---|
| 277 | ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock)) |
---|
| 278 | ! enddo |
---|
| 279 | ! enddo |
---|
[3178] | 280 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
| 281 | call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
[3149] | 282 | else |
---|
[3532] | 283 | ! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth |
---|
[3149] | 284 | tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope) |
---|
[3178] | 285 | call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) |
---|
| 286 | call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) |
---|
[3149] | 287 | tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope) |
---|
[3178] | 288 | call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope)) |
---|
[3203] | 289 | |
---|
[3532] | 290 | do iloop = nsoil_PCM + 1,nsoil_PEM |
---|
[3149] | 291 | tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope) |
---|
| 292 | enddo |
---|
| 293 | endif !found |
---|
[2794] | 294 | |
---|
[3149] | 295 | do it = 1,timelen |
---|
[3532] | 296 | tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) |
---|
[3149] | 297 | enddo |
---|
[3532] | 298 | do isoil = nsoil_PCM + 1,nsoil_PEM |
---|
| 299 | 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) |
---|
[3149] | 300 | enddo |
---|
| 301 | enddo ! islope |
---|
| 302 | write(*,*) 'PEMETAT0: SOIL TEMP done' |
---|
[2794] | 303 | |
---|
[3076] | 304 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2794] | 305 | !3. Ice Table |
---|
[3493] | 306 | if (icetable_equilibrium) then |
---|
| 307 | call get_field("ice_table_depth",ice_table_depth,found) |
---|
[3170] | 308 | if (.not. found) then |
---|
[3493] | 309 | write(*,*)'PEM settings: failed loading <ice_table_depth>' |
---|
[3170] | 310 | write(*,*)'will reconstruct the values of the ice table given the current state' |
---|
[3532] | 311 | call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) |
---|
[3525] | 312 | call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
[3170] | 313 | do islope = 1,nslope |
---|
[3178] | 314 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
[3170] | 315 | enddo |
---|
| 316 | endif |
---|
| 317 | write(*,*) 'PEMETAT0: ICE TABLE done' |
---|
[3493] | 318 | else if (icetable_dynamic) then |
---|
[3532] | 319 | call get_field("ice_porefilling",ice_porefilling,found) |
---|
| 320 | if (.not. found) then |
---|
| 321 | write(*,*)'PEM settings: failed loading <ice_porefilling>' |
---|
| 322 | ice_porefilling = 0. |
---|
| 323 | endif |
---|
[3493] | 324 | call get_field("ice_table_depth",ice_table_depth,found) |
---|
| 325 | if (.not. found) then |
---|
| 326 | write(*,*)'PEM settings: failed loading <ice_table_depth>' |
---|
| 327 | write(*,*)'will reconstruct the values of the ice table given the current state' |
---|
[3512] | 328 | ice_table_depth = -9999. |
---|
[3525] | 329 | call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
[3493] | 330 | do islope = 1,nslope |
---|
| 331 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
| 332 | enddo |
---|
| 333 | endif |
---|
| 334 | write(*,*) 'PEMETAT0: ICE TABLE done' |
---|
[3149] | 335 | endif |
---|
[2794] | 336 | |
---|
[3076] | 337 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2863] | 338 | !4. CO2 & H2O Adsorption |
---|
[3149] | 339 | if (adsorption_pem) then |
---|
| 340 | do islope = 1,nslope |
---|
| 341 | write(num,fmt='(i2.2)') islope |
---|
| 342 | call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found) |
---|
| 343 | if (.not. found) then |
---|
| 344 | m_co2_regolith_phys = 0. |
---|
| 345 | exit |
---|
| 346 | endif |
---|
| 347 | enddo |
---|
| 348 | do islope=1,nslope |
---|
| 349 | write(num,fmt='(i2.2)') islope |
---|
| 350 | call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2) |
---|
| 351 | if (.not. found2) then |
---|
| 352 | m_h2o_regolith_phys = 0. |
---|
| 353 | exit |
---|
| 354 | endif |
---|
| 355 | enddo |
---|
[2794] | 356 | |
---|
[3498] | 357 | call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, & |
---|
[3149] | 358 | m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys) |
---|
[3076] | 359 | |
---|
[3149] | 360 | if (.not. found) deltam_co2_regolith_phys = 0. |
---|
[3456] | 361 | if (.not. found2) deltam_h2o_regolith_phys = 0. |
---|
[2849] | 362 | |
---|
[3149] | 363 | write(*,*) 'PEMETAT0: CO2 & H2O adsorption done' |
---|
| 364 | endif ! adsorption_pem |
---|
| 365 | endif ! soil_pem |
---|
[3076] | 366 | |
---|
[3149] | 367 | call close_startphy |
---|
[2779] | 368 | |
---|
[3493] | 369 | else !No startpem, let's build all by hand |
---|
[2794] | 370 | |
---|
[3297] | 371 | write(*,*)'There is no "startpem.nc"!' |
---|
| 372 | |
---|
[3149] | 373 | ! h2o ice |
---|
| 374 | h2o_ice = 0. |
---|
[3214] | 375 | write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.' |
---|
| 376 | do ig = 1,ngrid |
---|
[3308] | 377 | if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice |
---|
[3214] | 378 | enddo |
---|
[2888] | 379 | |
---|
[3149] | 380 | ! co2 ice |
---|
[3214] | 381 | write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.' |
---|
[3149] | 382 | co2_ice = perennial_co2ice |
---|
[2794] | 383 | |
---|
[3430] | 384 | ! Stratification (layerings) |
---|
| 385 | nb_str_max = 1 |
---|
| 386 | if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.' |
---|
| 387 | |
---|
[3039] | 388 | if (soil_pem) then |
---|
[2794] | 389 | |
---|
| 390 | !a) Thermal inertia |
---|
[3149] | 391 | do islope = 1,nslope |
---|
| 392 | do ig = 1,ngrid |
---|
| 393 | if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then |
---|
| 394 | !!! transition |
---|
| 395 | delta = depth_breccia |
---|
| 396 | TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & |
---|
| 397 | (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ & |
---|
| 398 | ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) |
---|
| 399 | do iloop=index_breccia + 2,index_bedrock |
---|
| 400 | TI_PEM(ig,iloop,islope) = TI_breccia |
---|
| 401 | enddo |
---|
| 402 | else ! we keep the high ti values |
---|
| 403 | do iloop=index_breccia + 1,index_bedrock |
---|
| 404 | TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope) |
---|
| 405 | enddo |
---|
| 406 | endif |
---|
| 407 | !! transition |
---|
| 408 | delta = depth_bedrock |
---|
| 409 | TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ & |
---|
| 410 | (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ & |
---|
| 411 | ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2)))) |
---|
| 412 | do iloop = index_bedrock + 2,nsoil_PEM |
---|
| 413 | TI_PEM(ig,iloop,islope) = TI_bedrock |
---|
| 414 | enddo |
---|
| 415 | enddo |
---|
| 416 | enddo |
---|
[2794] | 417 | |
---|
[3149] | 418 | do iloop = 1,nsoil_PCM |
---|
| 419 | inertiedat_PEM(:,iloop) = inertiedat(:,iloop) |
---|
| 420 | enddo |
---|
| 421 | !!! zone de transition |
---|
| 422 | delta = depth_breccia |
---|
| 423 | do ig = 1,ngrid |
---|
| 424 | if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then |
---|
| 425 | inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & |
---|
| 426 | (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & |
---|
| 427 | ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2)))) |
---|
| 428 | do iloop = index_breccia + 2,index_bedrock |
---|
| 429 | inertiedat_PEM(ig,iloop) = TI_breccia |
---|
| 430 | enddo |
---|
| 431 | else |
---|
| 432 | do iloop = index_breccia + 1,index_bedrock |
---|
| 433 | inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia) |
---|
| 434 | enddo |
---|
| 435 | endif |
---|
| 436 | enddo |
---|
[2794] | 437 | |
---|
[3149] | 438 | !!! zone de transition |
---|
| 439 | delta = depth_bedrock |
---|
| 440 | do ig = 1,ngrid |
---|
| 441 | inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & |
---|
| 442 | (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + & |
---|
| 443 | ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2)))) |
---|
| 444 | enddo |
---|
[2779] | 445 | |
---|
[3149] | 446 | do iloop = index_bedrock + 2,nsoil_PEM |
---|
| 447 | do ig = 1,ngrid |
---|
| 448 | inertiedat_PEM(ig,iloop) = TI_bedrock |
---|
| 449 | enddo |
---|
| 450 | enddo |
---|
[2779] | 451 | |
---|
[3149] | 452 | write(*,*) 'PEMETAT0: TI done' |
---|
[2779] | 453 | |
---|
[3076] | 454 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 455 | !b) Soil temperature |
---|
[3149] | 456 | do islope = 1,nslope |
---|
| 457 | ! do ig = 1,ngrid |
---|
| 458 | ! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa |
---|
| 459 | ! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1)) |
---|
[2895] | 460 | ! |
---|
[3149] | 461 | ! do iloop=index_breccia+2,index_bedrock |
---|
| 462 | ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa |
---|
| 463 | ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia)) |
---|
| 464 | ! enddo |
---|
| 465 | ! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa |
---|
| 466 | ! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) |
---|
| 467 | ! |
---|
| 468 | ! do iloop=index_bedrock+2,nsoil_PEM |
---|
| 469 | ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa |
---|
| 470 | ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock)) |
---|
| 471 | ! enddo |
---|
| 472 | ! enddo |
---|
[3178] | 473 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
| 474 | call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
[2794] | 475 | |
---|
[3149] | 476 | ! First raw initialization |
---|
| 477 | do it = 1,timelen |
---|
| 478 | do isoil = nsoil_PCM + 1,nsoil_PEM |
---|
| 479 | tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope) |
---|
| 480 | enddo |
---|
| 481 | enddo |
---|
[3076] | 482 | |
---|
[3149] | 483 | do it = 1,timelen |
---|
| 484 | do isoil = nsoil_PCM + 1,nsoil_PEM |
---|
[3178] | 485 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it)) |
---|
[3149] | 486 | enddo |
---|
| 487 | enddo |
---|
[2794] | 488 | |
---|
[3149] | 489 | do isoil = nsoil_PCM + 1,nsoil_PEM |
---|
| 490 | do ig = 1,ngrid |
---|
| 491 | 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) |
---|
| 492 | enddo |
---|
| 493 | enddo |
---|
| 494 | enddo !islope |
---|
| 495 | write(*,*) 'PEMETAT0: TSOIL done' |
---|
[3039] | 496 | |
---|
[3076] | 497 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[3039] | 498 | !c) Ice table |
---|
[3493] | 499 | if (icetable_equilibrium) then |
---|
| 500 | call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) |
---|
[3525] | 501 | call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
[3170] | 502 | do islope = 1,nslope |
---|
[3178] | 503 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
[3170] | 504 | enddo |
---|
| 505 | write(*,*) 'PEMETAT0: Ice table done' |
---|
[3493] | 506 | else if (icetable_dynamic) then |
---|
[3532] | 507 | ice_porefilling = 0. |
---|
[3512] | 508 | ice_table_depth = -9999. |
---|
[3525] | 509 | call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
[3493] | 510 | do islope = 1,nslope |
---|
| 511 | call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) |
---|
| 512 | enddo |
---|
| 513 | write(*,*) 'PEMETAT0: Ice table done' |
---|
[3170] | 514 | endif |
---|
[3532] | 515 | |
---|
[3076] | 516 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[2794] | 517 | !d) Regolith adsorbed |
---|
[3149] | 518 | if (adsorption_pem) then |
---|
| 519 | m_co2_regolith_phys = 0. |
---|
| 520 | m_h2o_regolith_phys = 0. |
---|
[3498] | 521 | call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, & |
---|
[3149] | 522 | m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys) |
---|
| 523 | deltam_co2_regolith_phys = 0. |
---|
| 524 | deltam_h2o_regolith_phys = 0. |
---|
| 525 | endif |
---|
[2863] | 526 | |
---|
[3149] | 527 | write(*,*) 'PEMETAT0: CO2 adsorption done' |
---|
| 528 | endif !soil_pem |
---|
[2893] | 529 | endif ! of if (startphy_file) |
---|
| 530 | |
---|
[3076] | 531 | if (soil_pem) then ! Sanity check |
---|
| 532 | do ig = 1,ngrid |
---|
| 533 | do islope = 1,nslope |
---|
| 534 | do iloop = 1,nsoil_PEM |
---|
| 535 | if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1) |
---|
| 536 | enddo |
---|
| 537 | enddo |
---|
| 538 | enddo |
---|
| 539 | endif !soil_pem |
---|
[2794] | 540 | |
---|
[3076] | 541 | END SUBROUTINE pemetat0 |
---|
| 542 | |
---|
| 543 | END MODULE pemetat0_mod |
---|