Changeset 3031 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- Aug 24, 2023, 3:55:47 PM (18 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 3 edited
-
evol_h2o_ice_s_mod.F90 (modified) (5 diffs)
-
ice_table_mod.F90 (modified) (1 diff)
-
pem.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90
r3030 r3031 5 5 CONTAINS 6 6 7 SUBROUTINE evol_h2o_ice_s(qsurf,tendencies_h2o_ice_phys,& 8 iim_input,jjm_input,ngrid,cell_area,STOPPING,nslope) 7 SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING) 9 8 10 9 use temps_mod_evol, only: dt_pem … … 26 25 ! INPUT 27 26 28 INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope ! # of grid points along longitude/latitude/ total 29 REAL, intent(in) :: cell_area(ngrid) ! Area of each mesh grid 27 INTEGER, intent(in) :: ngrid ! # of grid points along longitude/latitude grid; 28 INTEGER, intent(in) :: nslope ! # of subslope 29 REAL, intent(in) :: cell_area(ngrid) ! Area of each mesh grid (m^2) 30 REAL, intent(in) :: delta_h2o_adsorbded(ngrid) ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2) 31 REAL, intent(in) :: delta_h2o_icetablesublim(ngrid) ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 30 32 31 33 ! OUTPUT 32 REAL, INTENT(INOUT) :: qsurf(ngrid,nslope) ! physical point field : Previous and actual density of water ice 34 REAL, INTENT(INOUT) :: qsurf(ngrid,nslope) ! physical point field : Previous and actual density of water ice (kg/m^2) 35 REAL, intent(inout) :: tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year (kg/m^2/year) 33 36 LOGICAL, INTENT(INOUT) :: STOPPING ! Stopping criterion 34 REAL, intent(inout) :: tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year35 37 36 38 ! local: 37 39 ! ---- 38 40 39 INTEGER :: i,j,i g0,islope! loop variable41 INTEGER :: i,j,islope ! loop variable 40 42 REAL :: pos_tend, neg_tend, real_coefficient,negative_part ! Variable to conserve water 41 43 REAL :: new_tendencies(ngrid,nslope) ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done … … 50 52 ! We compute the amount of water accumulating and sublimating 51 53 do i=1,ngrid 54 if(delta_h2o_adsorbded(i).GT.0) then 55 pos_tend=pos_tend+delta_h2o_adsorbded(i)*cell_area(i) 56 else 57 neg_tend=neg_tend+delta_h2o_adsorbded(i)*cell_area(i) 58 endif 59 if(delta_h2o_icetablesublim(i).GT.0) then 60 pos_tend=pos_tend+delta_h2o_icetablesublim(i)*cell_area(i) 61 else 62 neg_tend=neg_tend+delta_h2o_icetablesublim(i)*cell_area(i) 63 endif 52 64 do islope=1,nslope 53 65 if (qsurf(i,islope).GT.0) then … … 101 113 ! We compute the amount of water that is sublimated in excess 102 114 if (qsurf(i,islope).lt.0) then 103 negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope) 115 negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.) 104 116 qsurf(i,islope)=0. 105 117 tendencies_h2o_ice_phys(i,islope)=0. … … 119 131 do islope=1, nslope 120 132 if(new_tendencies(i,islope).GT.0) then ! In the place of accumulation of ice, we remove a bit of ice in order to conserve water 121 qsurf(i,islope)=qsurf(i,islope)-new_tendencies(i,islope)*real_coefficient*dt_pem 133 qsurf(i,islope)=qsurf(i,islope)-new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.) 122 134 endif 123 135 enddo -
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r2985 r3031 137 137 !!! -------------------------------------- 138 138 139 140 SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsoil,delta_m_h2o) 141 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 142 !!! 143 !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed 144 !!! 145 !!! Author: LL 146 !!! 147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 148 use comsoil_h_PEM, only: mlayer_PEM 149 use comslope_mod, only: subslope_dist,def_slope_mean 150 use comconst_mod,only: pi 151 implicit none 152 ! inputs 153 ! ----------- 154 integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 155 real,intent(in) :: former_ice_table_thickness(ngrid,nslope) ! ice table thickness at the former iteration [m] 156 real,intent(in) :: new_ice_table_thickness(ngrid,nslope) ! ice table thickness at the current iteration [m] 157 real,intent(out) :: ice_table_depth(ngrid,nslope) ! ice table depth [m] 158 real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope) ! Soil temperature [K] 159 ! outputs 160 ! ----------- 161 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] 162 163 ! local 164 real :: rho(ngrid,nslope) ! density of water ice [kg/m^3] 165 integer :: ig,islope,ilay,iref ! loop index 166 real :: Tice(ngrid,nslope) ! ice temperature [k] 167 ! Code 168 ! ----------- 169 rho(:,:) = 0. 170 Tice(:,:) = 0. 171 !1. First let's compute Tice using a linear interpolation between the mlayer level 172 do ig = 1,ngrid 173 do islope = 1,nslope 174 if(ice_table_depth(ig,islope).gt.0.) then 175 if (ice_table_depth(ig,islope).lt. 1e-10) then 176 Tice(ig,islope) = tsoil(ig,1,islope) 177 else 178 iref=0 ! initialize iref 179 do ilay=1,nsoil_PEM ! loop on layers to find the beginning of the ice table 180 if (ice_table_depth(ig,islope).ge.mlayer_PEM(ilay)) then 181 iref=ilay ! pure regolith layer up to here 182 else 183 ! correct iref was obtained in previous cycle 184 exit 185 endif 186 enddo 187 if(iref.eq.nsoil_PEM) then 188 Tice(ig,islope) = tsoil(ig,nsoil_PEM,islope) 189 else 190 Tice(ig,islope) = (tsoil(ig,iref+1,islope) - tsoil(ig,iref,islope))/(mlayer_PEM(iref+1) - mlayer_PEM(iref))* & 191 (ice_table_depth(ig,islope) - mlayer_PEM(iref)) + tsoil(ig,iref,islope) 192 endif 193 rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 194 endif 195 endif 196 enddo 197 enddo 198 199 200 !2. Let's compute the amount of ice that has sublimated in each subslope 201 do ig = 1,ngrid 202 do islope = 1,nslope 203 delta_m_h2o(ig) = delta_m_h2o(ig) + rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses 204 *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) 205 enddo 206 enddo 207 208 return 209 end subroutine 210 211 !!! -------------------------------------- 212 139 213 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) 140 214 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3030 r3031 57 57 use recomp_orb_param_mod, only: recomp_orb_param 58 58 use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, & 59 ini_ice_table_porefilling, computeice_table_equilibrium 59 ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi 60 60 use soil_thermalproperties_mod, only: update_soil_thermalproperties 61 61 … … 208 208 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 209 209 logical :: bool_sublim ! logical to check if there is sublimation or not 210 210 real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] 211 real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] 211 212 ! Some variables for the PEM run 212 213 real, parameter :: year_step = 1 ! timestep for the pem … … 499 500 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 500 501 allocate(delta_co2_adsorbded(ngrid)) 502 allocate(porefillingice_thickness_prev_iter(ngrid,nslope)) 503 allocate(delta_h2o_icetablesublim(ngrid)) 501 504 allocate(delta_h2o_adsorbded(ngrid)) 502 505 allocate(vmr_co2_pem_phys(ngrid,timelen)) … … 622 625 co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) 623 626 627 delta_h2o_icetablesublim(:) = 0. 628 624 629 do ig = 1,ngrid 625 630 do islope = 1,nslope … … 752 757 !------------------------ 753 758 write(*,*) "Evolution of h2o ice" 754 call evol_h2o_ice_s( qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)759 call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water) 755 760 756 761 write(*,*) "Evolution of co2 ice" … … 891 896 892 897 ! II_d.3 Update the ice table 898 porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:) 893 899 call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) 900 901 call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, & 902 delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere. 903 894 904 write(*,*) "Update soil propreties" 895 905 … … 1163 1173 deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved) 1164 1174 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) 1165 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys )1175 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) 1166 1176 !----------------------------- END OUTPUT ------------------------------ 1167 1177
Note: See TracChangeset
for help on using the changeset viewer.
