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