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