source: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90 @ 3296

Last change on this file since 3296 was 3206, checked in by jbclement, 10 months ago

PEM:
Cleanings of unused variables/arguments and bad type conversions.
JBC

File size: 6.5 KB
Line 
1MODULE evol_ice_mod
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
10
11use time_evol_mod, only: dt_pem
12
13implicit none
14
15!=======================================================================
16!
17! Routine to compute the evolution of CO2 ice
18!
19!=======================================================================
20!   arguments:
21!   ----------
22!   INPUT
23integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
24
25!   OUTPUT
26real, dimension(ngrid,nslope), intent(inout) :: co2_ice      ! Previous and actual density of CO2 ice
27real, 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
34write(*,*) 'Evolution of co2 ice'
35co2_ice = co2_ice + tend_co2_ice*dt_pem
36where (co2_ice < 0.)
37    co2_ice = 0.
38    tend_co2_ice = 0.
39end where
40
41END SUBROUTINE evol_co2_ice
42
43!=======================================================================
44
45SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
46
47use time_evol_mod,          only: dt_pem
48use 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
56implicit none
57
58!=======================================================================
59!
60! Routine that compute the evolution of the h2o ice
61!
62!=======================================================================
63!   arguments:
64!   ----------
65!   INPUT
66integer,                intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
67real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
68real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded      ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2)
69real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
70
71!   OUTPUT
72real, dimension(ngrid,nslope), intent(inout) :: h2o_ice      ! Previous and actual density of h2o ice (kg/m^2)
73real, dimension(ngrid,nslope), intent(inout) :: tend_h2o_ice ! Evolution of perennial ice over one year (kg/m^2/year)
74integer,                       intent(inout) :: stopPEM      ! Stopping criterion code
75
76!   local:
77!   ------
78integer                       :: i, islope                                           ! Loop variables
79real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
80real, dimension(ngrid,nslope) :: new_tendencies                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
81!=======================================================================
82if (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.)
153else ! ngrid == 1, i.e. in 1D
154    h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
155endif
156
157END SUBROUTINE evol_h2o_ice
158
159END MODULE evol_ice_mod
Note: See TracBrowser for help on using the repository browser.