| 1 | MODULE evol_ice_mod |
|---|
| 2 | |
|---|
| 3 | implicit none |
|---|
| 4 | |
|---|
| 5 | !======================================================================= |
|---|
| 6 | contains |
|---|
| 7 | !======================================================================= |
|---|
| 8 | |
|---|
| 9 | SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) |
|---|
| 10 | |
|---|
| 11 | use time_evol_mod, only: dt |
|---|
| 12 | use glaciers_mod, only: rho_co2ice |
|---|
| 13 | |
|---|
| 14 | implicit none |
|---|
| 15 | |
|---|
| 16 | !======================================================================= |
|---|
| 17 | ! |
|---|
| 18 | ! Routine to compute the evolution of CO2 ice |
|---|
| 19 | ! |
|---|
| 20 | !======================================================================= |
|---|
| 21 | ! arguments: |
|---|
| 22 | ! ---------- |
|---|
| 23 | ! INPUT |
|---|
| 24 | integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes |
|---|
| 25 | |
|---|
| 26 | ! OUTPUT |
|---|
| 27 | real, dimension(ngrid,nslope), intent(inout) :: co2_ice ! Previous and actual density of CO2 ice |
|---|
| 28 | real, dimension(ngrid,nslope), intent(inout) :: d_co2ice ! Evolution of perennial ice over one year |
|---|
| 29 | real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m] |
|---|
| 30 | |
|---|
| 31 | ! local: |
|---|
| 32 | ! ------ |
|---|
| 33 | real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice |
|---|
| 34 | !======================================================================= |
|---|
| 35 | ! Evolution of CO2 ice for each physical point |
|---|
| 36 | write(*,*) '> Evolution of CO2 ice' |
|---|
| 37 | |
|---|
| 38 | co2_ice_old = co2_ice |
|---|
| 39 | co2_ice = co2_ice + d_co2ice*dt |
|---|
| 40 | where (co2_ice < 0.) |
|---|
| 41 | co2_ice = 0. |
|---|
| 42 | d_co2ice = -co2_ice_old/dt |
|---|
| 43 | end where |
|---|
| 44 | |
|---|
| 45 | zshift_surf = d_co2ice*dt/rho_co2ice |
|---|
| 46 | |
|---|
| 47 | END SUBROUTINE evol_co2_ice |
|---|
| 48 | |
|---|
| 49 | !======================================================================= |
|---|
| 50 | |
|---|
| 51 | SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) |
|---|
| 52 | |
|---|
| 53 | use time_evol_mod, only: dt |
|---|
| 54 | use comslope_mod, only: subslope_dist, def_slope_mean |
|---|
| 55 | use glaciers_mod, only: rho_h2oice |
|---|
| 56 | |
|---|
| 57 | #ifndef CPP_STD |
|---|
| 58 | use comcstfi_h, only: pi |
|---|
| 59 | #else |
|---|
| 60 | use comcstfi_mod, only: pi |
|---|
| 61 | #endif |
|---|
| 62 | |
|---|
| 63 | implicit none |
|---|
| 64 | |
|---|
| 65 | !======================================================================= |
|---|
| 66 | ! |
|---|
| 67 | ! Routine to compute the evolution of h2o ice |
|---|
| 68 | ! |
|---|
| 69 | !======================================================================= |
|---|
| 70 | ! arguments: |
|---|
| 71 | ! ---------- |
|---|
| 72 | ! INPUT |
|---|
| 73 | integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes |
|---|
| 74 | real, dimension(ngrid), intent(in) :: cell_area ! Area of each mesh grid (m^2) |
|---|
| 75 | real, dimension(ngrid), intent(in) :: delta_h2o_adsorbed ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2) |
|---|
| 76 | real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) |
|---|
| 77 | |
|---|
| 78 | ! OUTPUT |
|---|
| 79 | real, dimension(ngrid,nslope), intent(inout) :: h2o_ice ! Previous and actual density of h2o ice (kg/m^2) |
|---|
| 80 | real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year) |
|---|
| 81 | integer, intent(inout) :: stopPEM ! Stopping criterion code |
|---|
| 82 | real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m] |
|---|
| 83 | |
|---|
| 84 | ! local: |
|---|
| 85 | ! ------ |
|---|
| 86 | integer :: i, islope ! Loop variables |
|---|
| 87 | real :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o |
|---|
| 88 | real, dimension(ngrid,nslope) :: new_tend ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done |
|---|
| 89 | !======================================================================= |
|---|
| 90 | write(*,*) '> Evolution of H2O ice' |
|---|
| 91 | if (ngrid /= 1) then ! Not in 1D |
|---|
| 92 | ! We compute the amount of condensing and sublimating h2o ice |
|---|
| 93 | pos_tend = 0. |
|---|
| 94 | neg_tend = 0. |
|---|
| 95 | do i = 1,ngrid |
|---|
| 96 | if (delta_h2o_adsorbed(i) > 0.) then |
|---|
| 97 | pos_tend = pos_tend + delta_h2o_adsorbed(i)*cell_area(i) |
|---|
| 98 | else |
|---|
| 99 | neg_tend = neg_tend + delta_h2o_adsorbed(i)*cell_area(i) |
|---|
| 100 | endif |
|---|
| 101 | if (delta_h2o_icetablesublim(i) > 0.) then |
|---|
| 102 | pos_tend = pos_tend + delta_h2o_icetablesublim(i)*cell_area(i) |
|---|
| 103 | else |
|---|
| 104 | neg_tend = neg_tend + delta_h2o_icetablesublim(i)*cell_area(i) |
|---|
| 105 | endif |
|---|
| 106 | do islope = 1,nslope |
|---|
| 107 | if (h2o_ice(i,islope) > 0.) then |
|---|
| 108 | if (d_h2oice(i,islope) > 0.) then |
|---|
| 109 | pos_tend = pos_tend + d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.) |
|---|
| 110 | else |
|---|
| 111 | neg_tend = neg_tend - d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.) |
|---|
| 112 | endif |
|---|
| 113 | endif |
|---|
| 114 | enddo |
|---|
| 115 | enddo |
|---|
| 116 | |
|---|
| 117 | new_tend = 0. |
|---|
| 118 | if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then |
|---|
| 119 | write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!" |
|---|
| 120 | write(*,*) "Tendencies on ice sublimating =", neg_tend |
|---|
| 121 | write(*,*) "Tendencies on ice increasing =", pos_tend |
|---|
| 122 | write(*,*) "This can be due to the absence of h2o ice in the PCM run!" |
|---|
| 123 | stopPEM = 2 |
|---|
| 124 | else |
|---|
| 125 | ! We adapt the tendencies to conserve h2o and do only exchange between grid points |
|---|
| 126 | if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation |
|---|
| 127 | where (d_h2oice < 0.) ! We lower the sublimating rate by a coefficient |
|---|
| 128 | new_tend = d_h2oice*pos_tend/neg_tend |
|---|
| 129 | elsewhere ! We don't change the condensing rate |
|---|
| 130 | new_tend = d_h2oice |
|---|
| 131 | end where |
|---|
| 132 | else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation |
|---|
| 133 | where (d_h2oice < 0.) ! We don't change the sublimating rate |
|---|
| 134 | new_tend = d_h2oice |
|---|
| 135 | elsewhere ! We lower the condensing rate by a coefficient |
|---|
| 136 | new_tend = d_h2oice*neg_tend/pos_tend |
|---|
| 137 | end where |
|---|
| 138 | endif |
|---|
| 139 | endif |
|---|
| 140 | |
|---|
| 141 | ! Evolution of the h2o ice for each physical point |
|---|
| 142 | h2o_ice = h2o_ice + new_tend*dt |
|---|
| 143 | |
|---|
| 144 | ! We compute the amount of h2o that is sublimated in excess |
|---|
| 145 | negative_part = 0. |
|---|
| 146 | do i = 1,ngrid |
|---|
| 147 | do islope = 1,nslope |
|---|
| 148 | if (h2o_ice(i,islope) < 0.) then |
|---|
| 149 | negative_part = negative_part - h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.) |
|---|
| 150 | h2o_ice(i,islope) = 0. |
|---|
| 151 | d_h2oice(i,islope) = 0. |
|---|
| 152 | endif |
|---|
| 153 | enddo |
|---|
| 154 | enddo |
|---|
| 155 | |
|---|
| 156 | ! We compute a coefficient by which we should remove the ice that has been added to places |
|---|
| 157 | ! even if this ice was contributing to an unphysical negative amount of ice at other places |
|---|
| 158 | if (abs(pos_tend) < 1.e-10) then |
|---|
| 159 | real_coefficient = 0. |
|---|
| 160 | else |
|---|
| 161 | real_coefficient = negative_part/pos_tend |
|---|
| 162 | endif |
|---|
| 163 | ! In the place of accumulation of ice, we remove a bit of ice in order to conserve h2o |
|---|
| 164 | do islope = 1,nslope |
|---|
| 165 | do i = 1,ngrid |
|---|
| 166 | if (new_tend(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tend(i,islope)*real_coefficient*dt*cos(def_slope_mean(islope)*pi/180.) |
|---|
| 167 | enddo |
|---|
| 168 | enddo |
|---|
| 169 | else ! ngrid == 1, i.e. in 1D |
|---|
| 170 | h2o_ice = h2o_ice + d_h2oice*dt |
|---|
| 171 | endif |
|---|
| 172 | |
|---|
| 173 | zshift_surf = d_h2oice*dt/rho_h2oice |
|---|
| 174 | |
|---|
| 175 | END SUBROUTINE evol_h2o_ice |
|---|
| 176 | |
|---|
| 177 | END MODULE evol_ice_mod |
|---|
