1 | !w ================================================================================================================================= |
---|
2 | ! MODULE : erosion |
---|
3 | ! |
---|
4 | ! CONTACT : Haicheng Zhang (sysuzhaicheng@163.com), Ronny Lauerwald |
---|
5 | ! |
---|
6 | ! LICENCE : IPSL (2006) |
---|
7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
8 | ! |
---|
9 | !>\BRIEF This module calculating the daily water-erosion rate. |
---|
10 | !! |
---|
11 | !!\n DESCRIPTION: None |
---|
12 | !! |
---|
13 | !! RECENT CHANGE(S): None |
---|
14 | !! |
---|
15 | !! REFERENCE(S) : |
---|
16 | !! |
---|
17 | !! SVN : |
---|
18 | !! $HeadURL$ |
---|
19 | !! $Date$ |
---|
20 | !! $Revision$ |
---|
21 | !! \n |
---|
22 | !_ ================================================================================================================================ |
---|
23 | !!========================================================================================= |
---|
24 | !!! Latest changes ------- by ZHC |
---|
25 | !!! Time = 20170809 |
---|
26 | !!! Change contents: |
---|
27 | !!! 1) The eroded soil carbon and belowground litter is calculated based on the carbon and litter contents in the top 7 soil layers, |
---|
28 | !!! rather than based on the carbon and litter contents in the top 8 soil layers in previous version. |
---|
29 | !!========================================================================================= |
---|
30 | MODULE erosion |
---|
31 | ! |
---|
32 | ! |
---|
33 | ! routines called : restput, restget |
---|
34 | ! |
---|
35 | USE time, ONLY : one_day, dt_sechiba |
---|
36 | USE ioipsl |
---|
37 | USE xios_orchidee |
---|
38 | USE ioipsl_para |
---|
39 | USE constantes |
---|
40 | USE constantes_soil |
---|
41 | USE constantes_var |
---|
42 | USE vertical_soil_var |
---|
43 | USE pft_parameters |
---|
44 | USE sechiba_io_p |
---|
45 | USE slowproc |
---|
46 | USE interpol_help |
---|
47 | USE grid |
---|
48 | USE mod_orchidee_para |
---|
49 | |
---|
50 | IMPLICIT NONE |
---|
51 | |
---|
52 | ! public routines : |
---|
53 | |
---|
54 | PRIVATE |
---|
55 | PUBLIC :: erosion_initialize, erosion_main, erosion_clear |
---|
56 | |
---|
57 | ! variables shared by all subroutines inside the erosion module : declaration and initialisation |
---|
58 | LOGICAL, SAVE :: l_first_erosion=.TRUE. !! Logical to initialize the routine. |
---|
59 | !! Initialisation has to be done one |
---|
60 | !! time (true/false) |
---|
61 | REAL(r_std),SAVE :: dt_erosion !! Erosion time step (s) |
---|
62 | REAL(r_std),SAVE :: time_counter_ero !! Time counter (s) |
---|
63 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: gridarea !! grid area (m^2) |
---|
64 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: cveg_f !! vegetation cover factor, unitless, 0-1 |
---|
65 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: kflitter !! vegetation cover factor, unitless, 0-1 |
---|
66 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: kfroot !! vegetation cover factor, unitless, 0-1 |
---|
67 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: meanrunoff !! Accumulated runoff(kg/m^{-2} = mm) |
---|
68 | !! @tex $(kg m^{-2} dt^{-1})$ @endtex |
---|
69 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: peakrunoff !! Peak runoff over one day (mm/s) |
---|
70 | !! @tex $(kg m^{-2} dt^{-1})$ @endtex |
---|
71 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: meanrundra !! Accumulated runoff+drainage(kg/m^{-2} = mm) |
---|
72 | !! @tex $(kg m^{-2} dt^{-1})$ @endtex |
---|
73 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: peakrundra !! Peak (runoff +drainage) over one day (mm/s) |
---|
74 | !! @tex $(kg m^{-2} dt^{-1})$ @endtex |
---|
75 | INTEGER(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: no_events !! number of rainfall events |
---|
76 | !! @tex $(kg m^{-2} dt^{-1})$ @endtex |
---|
77 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: eroc_active !! Active SOC export (gC/m^2/day) |
---|
78 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: eroc_slow !! Slow SOC export (gC/m^2/day) |
---|
79 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: eroc_passive !! Passive SOC export (gC/m^2/day) |
---|
80 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tot_7c !! Total Carbon concentration in the top 7 soil layers (gC m^{-2}) |
---|
81 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tot_7litbelow !! Total belowground litter in the top 7 soil layers (gC m^{-2}) |
---|
82 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: tot_7doc !! Total DOC in the top 7 soil layers (gC m^{-2}) |
---|
83 | !!********************* End veraibles for history files *********************** |
---|
84 | CONTAINS |
---|
85 | !! ====================================================================================================================== |
---|
86 | !! SUBROUTINE : erosion_initialize |
---|
87 | !! Initialization of the erosion module |
---|
88 | !! Read the reference soil erosion rate |
---|
89 | !! ====================================================================================================================== |
---|
90 | !! 1. Initialize |
---|
91 | SUBROUTINE erosion_initialize(kjit,nbpt,index,rest_id, lalo,neighbours,resolution,contfrac,std_totsed,gridarea,effgrid_perc) |
---|
92 | !! |
---|
93 | IMPLICIT NONE |
---|
94 | ! LOGICAL, SAVE :: l_first_erosion=.TRUE. !! Logical to initialize the routine. |
---|
95 | !! 0.1 Input variables |
---|
96 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
97 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
98 | INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) |
---|
99 | INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
100 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) |
---|
101 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point |
---|
102 | !! (unitless; 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
103 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) |
---|
104 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
105 | !! Output variables |
---|
106 | REAL(r_std), INTENT(out) :: std_totsed(nbpt) !! standard sediment discharge amount rate for each ORCHIDEE cell(ton) |
---|
107 | REAL(r_std), INTENT(out) :: gridarea(nbpt) !! Area of each grid-cell () |
---|
108 | REAL(r_std), INTENT(out) :: effgrid_perc(nbpt) !! percentage of the effective area in a ORCHIDEE pixel (0-1) |
---|
109 | !! Local variables |
---|
110 | INTEGER(i_std) :: ig |
---|
111 | !! 1.1 First call only |
---|
112 | dt_erosion = one_day |
---|
113 | WRITE(numout,*) 'Start the first time calling of erosion module', l_first_erosion |
---|
114 | IF (l_first_erosion) THEN |
---|
115 | WRITE(numout,*) 'This is the first time calling of erosion module' |
---|
116 | ALLOCATE(cveg_f(nbpt,nvm)) |
---|
117 | cveg_f = zero |
---|
118 | ALLOCATE(kflitter(nbpt,nvm)) |
---|
119 | kflitter = zero |
---|
120 | ALLOCATE(kfroot(nbpt,nvm)) |
---|
121 | kfroot = zero |
---|
122 | ALLOCATE(peakrunoff(nbpt)) |
---|
123 | peakrunoff = zero |
---|
124 | ALLOCATE(meanrunoff(nbpt)) |
---|
125 | meanrunoff = zero |
---|
126 | ALLOCATE(peakrundra(nbpt)) |
---|
127 | peakrundra = zero |
---|
128 | ALLOCATE(meanrundra(nbpt)) |
---|
129 | meanrundra = zero |
---|
130 | ALLOCATE(eroc_active(nbpt,nvm)) |
---|
131 | eroc_active = zero |
---|
132 | ALLOCATE(eroc_slow(nbpt,nvm)) |
---|
133 | eroc_slow = zero |
---|
134 | ALLOCATE(eroc_passive(nbpt,nvm)) |
---|
135 | eroc_passive = zero |
---|
136 | ALLOCATE(tot_7c(nbpt,ncarb,nvm)) |
---|
137 | tot_7c = zero |
---|
138 | ALLOCATE(tot_7litbelow(nbpt,nlitt,nvm)) |
---|
139 | tot_7litbelow = zero |
---|
140 | ALLOCATE(tot_7doc(nbpt,nvm,ndoc,npool)) |
---|
141 | tot_7doc = zero |
---|
142 | ALLOCATE(no_events(nbpt)) |
---|
143 | no_events = 0 |
---|
144 | ! Reading input data for calculating erosion (flacc, fldir, std_totsed) |
---|
145 | WRITE(numout,*) 'Start read erosion map' |
---|
146 | CALL erosionmap(nbpt,index,lalo,neighbours,resolution,contfrac,std_totsed,gridarea,effgrid_perc) |
---|
147 | ! Test reading input data |
---|
148 | !DO ig=1,nbpt |
---|
149 | ! WRITE(numout,*) 'igrid', ig,'lat',lalo(ig,1),'lon',lalo(ig,2) |
---|
150 | !ENDDO |
---|
151 | ! WRITE(numout,*) 'std_totsed', std_totsed(:) |
---|
152 | ! WRITE(numout,*) 'gridarea', gridarea(:) |
---|
153 | ! WRITE(numout,*) 'effgrid_perc', effgrid_perc(:) |
---|
154 | ! |
---|
155 | l_first_erosion = .FALSE. |
---|
156 | WRITE(numout,*) 'End first_erosion' |
---|
157 | ENDIF !l_first_erosion |
---|
158 | END SUBROUTINE erosion_initialize |
---|
159 | !! ====================================================================================================================================== |
---|
160 | !! SUBROUTINE : erosion_main |
---|
161 | !! ====================================================================================================================================== |
---|
162 | SUBROUTINE erosion_main(kjit,nbpt,rest_id,dt_sechiba,index, ldrestart_read, ldrestart_write, & |
---|
163 | & lalo,neighbours,resolution,temp_sol,bulkdens,textfrac,laydep,contfrac,runoff,drainage,veget, & |
---|
164 | & veget_max,hist_id,std_totsed,gridarea,effgrid_perc,rootmass,litter_above,litter_below, & |
---|
165 | & carbon, DOC, lignin_struc_above,lignin_struc_below,depth_deepsoil,seddep_rate, & |
---|
166 | & sed_deposition_d, poc_deposition_d,totrunoff_d,peakrunoff_d,totrundra_d,peakrundra_d, & |
---|
167 | & cf_veget,kf_litter,kf_root,erocap,erodep,eroc,tot_erocap,ave_erodep,tot_eroc, & |
---|
168 | & POC_EXP_agg,DOC_ERO_agg,SED_EXP_agg) |
---|
169 | IMPLICIT NONE |
---|
170 | !WRITE(numout,*) 'Erosion simulation0', fldir(1238726) |
---|
171 | !! 0.1 Input variables |
---|
172 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
173 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
174 | INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
175 | REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff @tex $(kg m^{-2} dt^{-1}, mm/dt)$ @endtex |
---|
176 | REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage @tex $(kg m^{-2} dt^{-1}, mm/dt)$ @endtex |
---|
177 | REAL(r_std), INTENT(in) :: dt_sechiba !! Time step of the model (s) |
---|
178 | INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) |
---|
179 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) |
---|
180 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point |
---|
181 | !! (unitless; 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
182 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) |
---|
183 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
184 | INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless) |
---|
185 | REAL(r_std),INTENT (in) :: veget(nbpt,nvm) !! Fraction of vegetation type (unitless;0-1) |
---|
186 | REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1) |
---|
187 | LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read (true/false) |
---|
188 | LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write (true/false) |
---|
189 | REAL(r_std), INTENT(in) :: temp_sol(nbpt) !! Surface temperature (K) |
---|
190 | REAL(r_std), INTENT (in) :: bulkdens(nbpt) !! bulk_density (kg m-3) zhc |
---|
191 | REAL(r_std), INTENT (in) :: textfrac(nbpt,nctext) !! Fraction of different soil particle sizes (0-1) |
---|
192 | REAL(r_std), INTENT(in) :: std_totsed(nbpt) !! standard sediment discharge amount rate for each ORCHIDEE cell(ton) |
---|
193 | REAL(r_std), INTENT(in) :: effgrid_perc(nbpt) !! percentage of the effective area in a ORCHIDEE pixel (0-1) |
---|
194 | REAL(r_std), INTENT(in) :: gridarea(nbpt) !! area of each grid cell (m^2) |
---|
195 | !! |
---|
196 | !!0.2 modified variables |
---|
197 | REAL(r_std), DIMENSION (nbpt,nvm,nparts,nelements), INTENT (inout) :: rootmass !! belowground biomass |
---|
198 | !! @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
199 | REAL(r_std), DIMENSION (nbpt,nlitt,nvm,nelements), INTENT (inout) :: litter_above !! Above ground metabolic and structural litter |
---|
200 | !! @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
201 | REAL(r_std), DIMENSION (nbpt,nlitt,nvm,nslmd,nelements), INTENT (inout) :: litter_below !! Below ground metabolic and structural litter |
---|
202 | !! per ground area @tex $(gC m^{-2})$ !!ORCHIDEE grid |
---|
203 | REAL(r_std), DIMENSION (nbpt,ncarb,nvm,nslmd), INTENT (inout) :: carbon !! Soil carbon pools per ground area: active, slow, or |
---|
204 | !! passive, @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
205 | REAL(r_std), DIMENSION(nbpt,nvm,nslmd,ndoc,npool,nelements), INTENT(inout) :: DOC !! Dissolved Organic Carbon in soil |
---|
206 | !! The unit is given by m^2 of |
---|
207 | !! ground @tex $(gC m{-2} of ground)$ @endtex |
---|
208 | REAL(r_std), DIMENSION(nbpt,nvm), INTENT(inout) :: lignin_struc_above !! Ratio Lignin content in structural litter, |
---|
209 | !! above ground, (0-1, unitless) |
---|
210 | REAL(r_std), DIMENSION(nbpt,nvm,nslmd), INTENT(inout) :: lignin_struc_below !! Ratio Lignin content in structural litter, |
---|
211 | !! below ground, (0-1, unitless) |
---|
212 | REAL(r_std), DIMENSION(nbpt,nctext), INTENT(inout) :: sed_deposition_d !! Daily deposition of sediment at floodplains |
---|
213 | !! @tex $(kg m^2 day^{-1})$ @endtex |
---|
214 | REAL(r_std), DIMENSION(nbpt,ncarb), INTENT(inout) :: poc_deposition_d !! Daily deposition of POC at floodplains |
---|
215 | !! @tex $(kg m^2 day^{-1})$ @endtex |
---|
216 | REAL(r_std), DIMENSION(nbpt,nvm), INTENT (inout) :: depth_deepsoil !! Depth of the soil layer deeper than 2 m (m) |
---|
217 | !! |
---|
218 | !! 0.3 output variables |
---|
219 | REAL(r_std), DIMENSION(nslm),INTENT (out) :: laydep !! The lower limit of each (11) soil layer (m) ::zhc |
---|
220 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT (out) :: seddep_rate !! Sediment deposition rate (m day-1) for each PFT (m/day) |
---|
221 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: totrunoff_d !! daily total runoff (mm) |
---|
222 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: peakrunoff_d !! daily peak runoff rate.(mm/s = kg/m-2/s) |
---|
223 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: totrundra_d !! daily total runoff+drainage (mm) |
---|
224 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: peakrundra_d !! daily peak (runoff+drainage) rate.(mm/s = kg/m-2/s) |
---|
225 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT(out) :: erocap !! Erosion capacity from subgrid.(ton) |
---|
226 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT(out) :: erodep !! Erosion depth of sub_grid (m) under different PFTs (m) |
---|
227 | REAL(r_std), DIMENSION(nbpt,ncarb,nvm),INTENT(out) :: eroc !! Carbon export (gC/m^2/day) |
---|
228 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: tot_erocap !! Erosion capacity from subgrid.(ton) |
---|
229 | REAL(r_std), DIMENSION(nbpt),INTENT(out) :: ave_erodep !! Erosion depth of sub_grid (m) |
---|
230 | REAL(r_std), DIMENSION(nbpt,ncarb),INTENT(out) :: tot_eroc !! Carbon export (gC/m^2/day) |
---|
231 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT(out) :: cf_veget !! Daily vegetation cover factor |
---|
232 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT(out) :: kf_litter !! Daily aboveground litter cover factor |
---|
233 | REAL(r_std), DIMENSION(nbpt,nvm),INTENT(out) :: kf_root !! Daily belowground litter cover factor |
---|
234 | REAL(r_std), DIMENSION(nbpt,ncarb),INTENT(out) :: POC_EXP_agg !! POC exports along with surface runoff, in |
---|
235 | !! @tex $(gC m^{-2} day^{-1})$ @endtex |
---|
236 | REAL(r_std), DIMENSION(nbpt,nctext),INTENT(out) :: SED_EXP_agg !! sediment exports along with surface runoff, in |
---|
237 | !! @tex $(g m^{-2} day^{-1})$ @endtex |
---|
238 | REAL(r_std), DIMENSION(nbpt,ncarb),INTENT(out) :: DOC_ERO_agg !! DOC exports along with surface runoff, in |
---|
239 | !! @tex $(g m^{-2} day^{-1})$ @endtex |
---|
240 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
241 | ! |
---|
242 | !! 0.4 Local variables |
---|
243 | INTEGER(i_std) :: ig,jg,lg,fg,ng,ii,icarb,ilitt,ibdl !! integers for loops |
---|
244 | INTEGER(i_std) :: ilat,ilon,inbpt !! temperal variable for latitude and longitude |
---|
245 | INTEGER(i_std) :: max_flacc !! maximum value of flow accumulation |
---|
246 | REAL(r_std) :: erocap_jg !! temperal variable for erosion capacity |
---|
247 | REAL(r_std) :: meanrun_up,peakrun_up,meanrundra_up,peakrundra_up,cveg_up,kflitter_up,kfroot_up !! temperal variable for erosion capacity |
---|
248 | REAL(r_std) :: littotal,belowbiom !! vegetation cover fraction (0-1) |
---|
249 | REAL(r_std) :: vegf !! vegetation cover fraction (0-1) |
---|
250 | REAL(r_std) :: max_depth !! vegetation cover fraction (0-1) |
---|
251 | REAL(r_std) :: lignin_ave !! Average ratio of Lignin content in belowground structural litter in top 7 soil layers |
---|
252 | ! REAL(r_std) :: ave_c7,ave_lit !! temperal variable for erosion capacity |
---|
253 | ! REAL(r_std) :: nsubcell !! potentially maximum number of finer cells(erosion map) |
---|
254 | CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) |
---|
255 | CHARACTER(LEN=10) :: carb_suff !! To assign a suffix to the variables name indicating carbon type |
---|
256 | !! (active, slow, passive) |
---|
257 | CHARACTER(LEN=10) :: litt_suff !! To assign a suffix to the variables name indicating litterfall type |
---|
258 | !! (structual, metabolic) |
---|
259 | CHARACTER(LEN=10) :: solly_suff !! To store variables names for I/O soil layer |
---|
260 | REAL(r_std), DIMENSION(nbpt,ncarb,nvm,nslmd) ::carbon_old !! Soil carbon pools per ground area: active, slow, or |
---|
261 | !! passive, @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
262 | REAL(r_std), DIMENSION(nbpt,nlitt,nvm,nslmd) ::litbe_old !! belowground litter pools per ground area: |
---|
263 | !! @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
264 | REAL(r_std), DIMENSION(nbpt,nvm,nslmd) ::litbe_tot !! total belowground litter pools per ground area: |
---|
265 | !! @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
266 | REAL(r_std), DIMENSION(nbpt,nvm,nslmd,ndoc,npool) ::doc_old !! DOC pools per ground area: |
---|
267 | !! @tex $(gC m^{-2})$ @endtex !!ORCHIDEE grid |
---|
268 | REAL(r_std), DIMENSION(nbpt,nvm) ::erorate !! Daily average net soil loss rate (g m^{-2} dt_slow^{-1} PFT-1) |
---|
269 | REAL(r_std), DIMENSION(nbpt) ::sed_deposition_tot !! Daily total sediment deposition rate (g m^{-2} day-1) |
---|
270 | REAL(r_std), DIMENSION(nbpt,ncarb,nvm) ::soc_deposition !! Daily POC deposition rate (g C m^{-2} day-1 PFT-1) |
---|
271 | REAL(r_std), DIMENSION(nbpt) :: bio_frac !! total wet deposition of DOC onto the flooded |
---|
272 | REAL(r_std), DIMENSION(nslm) :: dly !! Depth of each soil layer (m) |
---|
273 | REAL(r_std), DIMENSION(nbpt,nlitt,nvm) :: erolitbe !! Eroded belowground litter (gC/m^2/day) |
---|
274 | REAL(r_std), DIMENSION(nbpt,nvm,ndoc,npool) :: erodoc !! Eroded doc (gC/m^2/day) |
---|
275 | REAL(r_std), DIMENSION(nbpt,nctext) :: sed_text !! Texture of sediment exported to the river (fast resvoir) |
---|
276 | !_ ================================================================ |
---|
277 | |
---|
278 | ! |
---|
279 | laydep(nslm) = maxsoildep |
---|
280 | DO ig = 1,nslm-1 |
---|
281 | laydep(ig) = maxsoildep/(2**(nslm-1)-1)*((2**(ig-1) -1)+(2**(ig)-1))/2.0 |
---|
282 | ENDDO |
---|
283 | |
---|
284 | dly(1)=laydep(1) |
---|
285 | DO ig = 2,nslm |
---|
286 | dly(ig)=laydep(ig)-laydep(ig-1) |
---|
287 | ENDDO |
---|
288 | ! WRITE(numout,*) 'Depth_Soillayer', laydep(:) |
---|
289 | |
---|
290 | !! 1.1 Aggregate runoff values per day. Extract the daily total runoff amount and peak runoff rate |
---|
291 | !! peak runoff rate = maximum of the 48 half-hour runoff rates winthin a day |
---|
292 | DO ig=1,nbpt |
---|
293 | ! Only runoff considered |
---|
294 | meanrunoff(ig) = meanrunoff(ig) + runoff(ig) |
---|
295 | IF (runoff(ig) .GT. peakrunoff(ig)) THEN |
---|
296 | peakrunoff(ig) = runoff(ig) |
---|
297 | ENDIF !(runoff(ig) .GT. peakrunoff(ig)) |
---|
298 | |
---|
299 | IF (runoff(ig) .GT. min_sechiba) THEN |
---|
300 | no_events(ig) = no_events(ig) + 1 |
---|
301 | ENDIF !(runoff(ig) .GT. min_sechiba) |
---|
302 | ! |
---|
303 | ! Both runoff and drainage considered |
---|
304 | meanrundra(ig) = meanrundra(ig) + runoff(ig)+drainage(ig) |
---|
305 | IF (runoff(ig)+drainage(ig) .GT. peakrundra(ig)) THEN |
---|
306 | peakrundra(ig) = runoff(ig)+drainage(ig) |
---|
307 | ENDIF !(runoff(ig) .GT. peakrunoff(ig)) |
---|
308 | ENDDO !ig=1,nbpt |
---|
309 | ! |
---|
310 | time_counter_ero = time_counter_ero + dt_sechiba |
---|
311 | ! WRITE(numout,*) 'time_counter_ero',time_counter_ero,'dt_sechiba',dt_sechiba |
---|
312 | ! WRITE(numout,*) 'runoff',runoff(1),'meanrunoff',meanrunoff(1),'peakrunoff',peakrunoff(1) |
---|
313 | ! |
---|
314 | seddep_rate = zero !! Sediment deposition rate (m day-1) for each PFT (m/day) |
---|
315 | totrunoff_d = zero !! daily total runoff (mm) |
---|
316 | peakrunoff_d = zero !! daily peak runoff rate.(mm/s = kg/m-2/s) |
---|
317 | totrundra_d = zero !! daily total runoff+drainage (mm) |
---|
318 | peakrundra_d = zero !! daily peak (runoff+drainage) rate.(mm/s = kg/m-2/s) |
---|
319 | erocap = zero !! Erosion capacity from subgrid.(ton) |
---|
320 | erodep= zero !! Erosion depth of sub_grid (m) under different PFTs (m) |
---|
321 | eroc = zero !! Carbon export (gC/m^2/day) |
---|
322 | tot_erocap = zero !! Erosion capacity from subgrid.(ton) |
---|
323 | ave_erodep = zero !! Erosion depth of sub_grid (m) |
---|
324 | tot_eroc = zero !! Carbon export (gC/m^2/day) |
---|
325 | cf_veget = zero !! Daily vegetation cover factor |
---|
326 | kf_litter = zero !! Daily aboveground litter cover factor |
---|
327 | kf_root = zero !! Daily belowground litter cover factor |
---|
328 | POC_EXP_agg = zero !! POC exports along with surface runoff, in ex $(gC m^{-2} dt_slow^{-1})$ @endtex |
---|
329 | SED_EXP_agg = zero !! sediment exports along with surface runoff, in @tex $(g m^{-2} dt_slow^{-1})$ @endtex |
---|
330 | DOC_ERO_agg = zero !! DOC exports along with surface runoff, in @tex $(g m^{-2} dt_slow^{-1})$ @endtex |
---|
331 | !! 2 If the simulation for the whole 48 half-hour intervals of a day have been finished, |
---|
332 | !! we calculate the erosion. |
---|
333 | !! 2.1 output daily mean runoff and the peak runoff rate |
---|
334 | ! |
---|
335 | IF ( NINT(time_counter_ero) .GE. NINT(dt_erosion) ) THEN |
---|
336 | |
---|
337 | bio_frac(:)=SUM(veget_max(:,2:13),DIM=2) |
---|
338 | |
---|
339 | peakrunoff(:) = peakrunoff(:) / dt_sechiba |
---|
340 | peakrundra(:) = peakrundra(:) / dt_sechiba |
---|
341 | !! peakrunoff(:) = peakrunoff(:) |
---|
342 | WHERE (peakrunoff(:) .LT. min_sechiba) |
---|
343 | peakrunoff(:) = zero |
---|
344 | ENDWHERE |
---|
345 | WHERE (peakrundra(:) .LT. min_sechiba) |
---|
346 | peakrundra(:) = zero |
---|
347 | ENDWHERE |
---|
348 | ! |
---|
349 | WHERE (meanrunoff(:) .LT. min_sechiba) |
---|
350 | meanrunoff(:) = zero |
---|
351 | ENDWHERE |
---|
352 | WHERE (meanrundra(:) .LT. min_sechiba) |
---|
353 | meanrundra(:) = zero |
---|
354 | ENDWHERE |
---|
355 | totrunoff_d(:)=meanrunoff(:) |
---|
356 | peakrunoff_d(:)=peakrunoff(:) |
---|
357 | totrundra_d(:)=meanrundra(:) |
---|
358 | peakrundra_d(:)=peakrundra(:) |
---|
359 | ! |
---|
360 | ! WRITE(numout,*) 'no of events: ', no_events(:) |
---|
361 | ! WRITE(numout,*) 'peak runoff was: ', peakrunoff(:) |
---|
362 | ! WRITE(numout,*) 'mean runoff was: ', meanrunoff(:) |
---|
363 | ! WRITE(numout,*) 'peak runoff was: ', peakrundra(:) |
---|
364 | ! WRITE(numout,*) 'mean runoff was: ', meanrundra(:) |
---|
365 | ! WRITE(numout,*) 'Erosion_bulkdens: ', bulkdens(:) |
---|
366 | ! WRITE(numout,*) 'Erosion_ZHC1' |
---|
367 | ! WRITE(numout,*) 'litt_ab: ', MAXVAL(litter_above),MINVAL(litter_above) |
---|
368 | ! WRITE(numout,*) 'litt_be: ', MAXVAL(litter_below),MINVAL(litter_below) |
---|
369 | ! WRITE(numout,*) 'ligninbe',MAXVAL(lignin_struc_below),MINVAL(lignin_struc_below) |
---|
370 | ! WRITE(numout,*) 'DOC',MAXVAL(DOC),MINVAL(DOC) |
---|
371 | ! WRITE(numout,*) 'SOC',MAXVAL(carbon),MINVAL(carbon) |
---|
372 | ! WRITE(numout,*) 'layer_depth1: ', laydep(7) |
---|
373 | ! |
---|
374 | ! 2.2 Calculate vegetation, litter, root factor |
---|
375 | ! 2.2.1 Calculate vegetation and management factor |
---|
376 | cf_veget(:,:)=0.0 |
---|
377 | cveg_f(:,:)=0.0 |
---|
378 | cf_veget(:,1)=zero |
---|
379 | cveg_f(:,1)=zero |
---|
380 | |
---|
381 | DO jg=1,nbpt |
---|
382 | cveg_f(jg,1)=0.0 |
---|
383 | DO lg=2,nvm |
---|
384 | IF (veget_max(jg,lg).GT.zero) THEN |
---|
385 | vegf=100.0*veget(jg,lg)/veget_max(jg,lg) |
---|
386 | IF (vegf.LT. 0.11) THEN |
---|
387 | cveg_f(jg,lg) = 1.0 |
---|
388 | ELSEIF (vegf .GT. 95.) THEN |
---|
389 | cveg_f(jg,lg) = min_cvegf |
---|
390 | ELSE |
---|
391 | cveg_f(jg,lg) = cvegf_coef1-cvegf_coef2*log10(vegf) |
---|
392 | ENDIF !(veget(jg,2:nvm) .LT. 0.11) |
---|
393 | ENDIF |
---|
394 | ENDDO !DO lg=2,nvm |
---|
395 | ENDDO !jg=1,nbpt |
---|
396 | cf_veget(:,:)=cveg_f(:,:) |
---|
397 | ! 2.2.2 Calculate factors of aboveground litter and belowground litter |
---|
398 | DO jg=1,nbpt |
---|
399 | kflitter(jg,1)=0.0 |
---|
400 | kfroot(jg,1)=0.0 |
---|
401 | DO lg=2,nvm |
---|
402 | IF (lg.LT.12) THEN |
---|
403 | littotal=0.001*(sum(litter_above(jg,:,lg,icarbon))+sum(litter_below(jg,:,lg,1:7,icarbon))) |
---|
404 | ELSE |
---|
405 | littotal=0.001*sum(litter_below(jg,:,lg,1:7,icarbon)) |
---|
406 | ENDIF |
---|
407 | ! |
---|
408 | IF (veget_max(jg,lg).GT.zero) THEN |
---|
409 | kflitter(jg,lg)=EXP(-littf_exponent*littotal) |
---|
410 | ELSE |
---|
411 | kflitter(jg,lg)=0.0 |
---|
412 | ENDIF |
---|
413 | ! |
---|
414 | belowbiom=0.001*(rootmass(jg,lg,isapbelow,icarbon)+ & |
---|
415 | & rootmass(jg,lg,iheartbelow,icarbon)+rootmass(jg,lg,iroot,icarbon)) |
---|
416 | IF (veget_max(jg,lg).GT.zero) THEN |
---|
417 | kfroot(jg,lg)=EXP(-rootf_exponent*belowbiom) |
---|
418 | ELSE |
---|
419 | kfroot(jg,lg)=0.0 |
---|
420 | ENDIF |
---|
421 | !WRITE(numout,*) 'ZHC2', jg,lg,littotal,kflitter(jg,lg) |
---|
422 | ENDDO !lg=2,nvm |
---|
423 | ENDDO !jg=1,nbpt |
---|
424 | kf_litter(:,:)=kflitter(:,:) |
---|
425 | kf_root(:,:)=kfroot(:,:) |
---|
426 | ! |
---|
427 | !! 2.3 Calculate erosion rate based on MUSLE equation |
---|
428 | ! 2.3.1 Calculate the sediment amount(ton) for each finer(sub-) grid (standard erosion grid) |
---|
429 | erocap(:,:)=zero |
---|
430 | tot_erocap(:)=zero |
---|
431 | erodep(:,:)=zero |
---|
432 | ave_erodep(:)=zero |
---|
433 | erorate(:,:)=zero |
---|
434 | |
---|
435 | DO ig=1,nbpt |
---|
436 | meanrun_up=totrunoff_d(ig)! |
---|
437 | peakrun_up=peakrunoff_d(ig)! |
---|
438 | IF (temp_sol(ig).LT.(273.15-5.0)) THEN ! H. Zhang: When soil temperature is lower than -5-degree, we assumed no erosion occure |
---|
439 | erocap(ig,:)=zero |
---|
440 | erodep(ig,:)=zero |
---|
441 | tot_erocap(ig)=zero |
---|
442 | ave_erodep(ig)=zero |
---|
443 | ELSE |
---|
444 | DO lg=2,nvm |
---|
445 | IF (veget_max(ig,lg).GT.zero) THEN |
---|
446 | ! Only surface runoff considered |
---|
447 | cveg_up=cf_veget(ig,lg) |
---|
448 | kflitter_up=kf_litter(ig,lg) |
---|
449 | kfroot_up=kf_root(ig,lg) |
---|
450 | erocap(ig,lg)=(meanrun_up*peakrun_up/(std_runoff*std_qp)/1000.)**runf_power & |
---|
451 | & *(cveg_up/std_cf)*std_totsed(ig)*kflitter_up*kfroot_up*veget_max(ig,lg) !*kfroot_up % unit ton/day (Mg/day) |
---|
452 | ! surface runoff+drainage |
---|
453 | ! erocap(ig)=(meanrundra_up*peakrundra_up/(std_runoff*std_qp)/1000)**0.56 & |
---|
454 | ! & *(cveg_up/std_cf)*std_totsed(ig)*kflitter_up !*kfroot_up |
---|
455 | IF (erocap(ig,lg).LT.zero) THEN |
---|
456 | erocap(ig,lg)=zero |
---|
457 | ENDIF |
---|
458 | ! depth of the sediment if all input sediments from upstream are deposited at one grid-cell |
---|
459 | ! Although soil erosion only occurs in the effective part. the erosion rate/depth is averaged to the whole grid cell |
---|
460 | IF (gridarea(ig).GT.zero.AND.bulkdens(ig).GT.zero) THEN |
---|
461 | erorate(ig,lg)=1.0E+6*erocap(ig,lg)/(gridarea(ig)*veget_max(ig,lg)) ! g soil /m^2 /day ##*effgrid_perc(ig) |
---|
462 | erodep(ig,lg)=1000.*erocap(ig,lg)/(gridarea(ig)*veget_max(ig,lg))/bulkdens(ig) !m/day ## *effgrid_perc(ig) |
---|
463 | !WRITE(numout,*) 'ZHC3_ero', ig,lg,effgrid_perc(ig),bulkdens(ig),erodep(ig,lg),erorate(ig,lg),std_totsed(ig) |
---|
464 | ENDIF |
---|
465 | ! |
---|
466 | ENDIF !(veget_max(ig,lg).GT.zero) THEN |
---|
467 | ENDDO !DO lg=1,nvm |
---|
468 | ! Average erosion rate for each grid-cell |
---|
469 | tot_erocap(ig)=SUM(erocap(ig,:)) ! ton/day |
---|
470 | IF (gridarea(ig).GT.zero.AND.bulkdens(ig).GT.zero) THEN |
---|
471 | ave_erodep(ig)=1000.*tot_erocap(ig)/(bio_frac(ig)*gridarea(ig))/bulkdens(ig) ! m/day ## *effgrid_perc(ig) |
---|
472 | ENDIF |
---|
473 | ENDIF !IF (temp_sol(ig).LT.(-2.0)) THEN |
---|
474 | ENDDO !DO lg=1,nbpt |
---|
475 | !WRITE(numout,*) 'erosion depth', SUM(ave_erodep(:)) |
---|
476 | ! |
---|
477 | !2.4 Calculate the processes of erosion rate of SOC |
---|
478 | ! Calculate from the upstream to outlets of the watershed based on flow accumulation |
---|
479 | !WRITE(numout,*) 'erosion_carbon', carbon(:,1,:,1) |
---|
480 | tot_eroc(:,:)=zero |
---|
481 | ! |
---|
482 | carbon_old(:,:,:,:)=carbon(:,:,:,:) |
---|
483 | litbe_old(:,:,:,:)=litter_below(:,:,:,:,icarbon) |
---|
484 | litbe_tot(:,:,:)=SUM(litbe_old(:,:,:,:),DIM=2) |
---|
485 | doc_old(:,:,:,:,:)=doc(:,:,:,:,:,icarbon) |
---|
486 | ! |
---|
487 | eroc(:,:,:)=zero |
---|
488 | erolitbe(:,:,:)=zero |
---|
489 | erodoc(:,:,:,:)=zero |
---|
490 | ! |
---|
491 | tot_7c(:,:,:)=zero |
---|
492 | tot_7litbelow(:,:,:)=zero |
---|
493 | tot_7doc(:,:,:,:)=zero |
---|
494 | |
---|
495 | !WRITE(numout,*) 'SUMero1_litter',SUM(litter_above),SUM(litter_below) |
---|
496 | !WRITE(numout,*) 'SUMero1_SOC',SUM(carbon) |
---|
497 | !WRITE(numout,*) 'SUMero1_DOC',SUM(DOC) |
---|
498 | !WRITE(numout,*) 'ligninbe',SUM(lignin_struc_below) |
---|
499 | !WRITE(numout,*) 'laydep',laydep(:) |
---|
500 | !WRITE(numout,*) 'dly',dly(:) |
---|
501 | |
---|
502 | DO ig=1,nbpt |
---|
503 | !WRITE(numout,*) 'SUMero1',ig,SUM(litter_below(ig,:,1,:,icarbon)),SUM(carbon(ig,:,1,:)), & |
---|
504 | ! lignin_struc_below(ig,1,:),SUM(DOC(ig,1,:,:,:,icarbon)) |
---|
505 | DO lg=2,nvm |
---|
506 | !WRITE(numout,*) 'SUMero2',ig,lg,SUM(litter_below(ig,:,lg,:,icarbon)),SUM(DOC(ig,1,:,:,:,icarbon)) |
---|
507 | |
---|
508 | erodep(ig,lg)=MIN(erodep(ig,lg),laydep(7)) ! H. Zhang: depth of daily eroded soil will not exceed laydep(7) (20 cm) |
---|
509 | |
---|
510 | IF (ISNAN(SUM(carbon(ig,3,lg,:)))) THEN |
---|
511 | carbon(ig,:,lg,:)=zero |
---|
512 | carbon_old(ig,:,lg,:)=zero |
---|
513 | ENDIF ! |
---|
514 | ! |
---|
515 | IF ((veget_max(ig,lg) .GT. 0).AND.(erodep(ig,lg).GT.0).AND.(MinVal(SUM(carbon(ig,:,lg,1:7),DIM=2)).GT.0)) THEN !!! 1<=fldir<=128 |
---|
516 | DO jg=1,7 |
---|
517 | tot_7c(ig,:,lg)=tot_7c(ig,:,lg)+carbon(ig,:,lg,jg) |
---|
518 | tot_7litbelow(ig,:,lg)=tot_7litbelow(ig,:,lg)+litter_below(ig,:,lg,jg,icarbon) |
---|
519 | ! tot_7doc(ig,lg,:,:)=tot_7doc(ig,lg,:,:)+DOC(ig,lg,jg,:,:,icarbon) |
---|
520 | lignin_ave= lignin_ave+lignin_struc_below(ig,lg,jg)*SUM(litter_below(ig,:,lg,jg,icarbon)) |
---|
521 | ENDDO !DO jg=,7 |
---|
522 | eroc(ig,:,lg) = erodep(ig,lg)/laydep(7)*tot_7c(ig,:,lg) |
---|
523 | erolitbe(ig,:,lg) = erodep(ig,lg)/laydep(7)*tot_7litbelow(ig,:,lg) |
---|
524 | !erodoc(ig,lg,:,:) = erodep(ig,lg)/laydep(7)*tot_7doc(ig,lg,:,:) |
---|
525 | |
---|
526 | tot_eroc(ig,:)=tot_eroc(ig,:)+ eroc(ig,:,lg)*veget_max(ig,lg) |
---|
527 | |
---|
528 | !! H. Zhang. Here we assumed that the soil erosion only takes sediment, POC and DOC into rivers. |
---|
529 | !! Belowground litter in the eroded soil layer will enter into the surface litter pools |
---|
530 | IF (SUM(tot_7litbelow(ig,:,lg)).GT.0) THEN |
---|
531 | lignin_ave=lignin_ave/SUM(tot_7litbelow(ig,:,lg)) |
---|
532 | ! lignin_struc_above(ig,lg)= (lignin_struc_above(ig,lg)*SUM(litter_above(ig,:,lg,icarbon))+ & |
---|
533 | ! lignin_ave*SUM(erolitbe(ig,:,lg)))/(SUM(litter_above(ig,:,lg,icarbon))+SUM(erolitbe(ig,:,lg))) |
---|
534 | litter_above(ig,:,lg,icarbon)=litter_above(ig,:,lg,icarbon) + erolitbe(ig,:,lg) |
---|
535 | ENDIF |
---|
536 | |
---|
537 | DO jg=1,7 ! top 7 soil layers |
---|
538 | WHERE (tot_7c(ig,:,lg).GT.zero) |
---|
539 | carbon(ig,:,lg,jg)=carbon_old(ig,:,lg,jg)-eroc(ig,:,lg)*carbon_old(ig,:,lg,jg)/tot_7c(ig,:,lg)+ & |
---|
540 | & (erodep(ig,lg)/dly(8)*carbon_old(ig,:,lg,8))*carbon_old(ig,:,lg,jg)/tot_7c(ig,:,lg) |
---|
541 | ENDWHERE |
---|
542 | |
---|
543 | IF (SUM(tot_7litbelow(ig,:,lg)).GT.zero) THEN |
---|
544 | WHERE (tot_7litbelow(ig,:,lg).GT.zero) |
---|
545 | litter_below(ig,:,lg,jg,icarbon)=litbe_old(ig,:,lg,jg)-erolitbe(ig,:,lg)*litbe_old(ig,:,lg,jg)/tot_7litbelow(ig,:,lg)+ & |
---|
546 | & (erodep(ig,lg)/dly(8)*litbe_old(ig,:,lg,8))*litbe_old(ig,:,lg,jg)/tot_7litbelow(ig,:,lg) |
---|
547 | ENDWHERE |
---|
548 | ! lignin_struc_below(ig,lg,jg)= (lignin_struc_below(ig,lg,jg)*(litbe_tot(ig,lg,jg)- & |
---|
549 | ! & SUM(erolitbe(ig,:,lg))*litbe_tot(ig,lg,jg)/SUM(tot_7litbelow(ig,:,lg)))+ & |
---|
550 | ! & lignin_struc_below(ig,lg,8)*(erodep(ig,lg)/dly(8)*litbe_tot(ig,lg,8))*litbe_tot(ig,lg,jg)/SUM(tot_7litbelow(ig,:,lg))) / & |
---|
551 | ! & SUM(litter_below(ig,:,lg,jg,icarbon)) |
---|
552 | ENDIF |
---|
553 | |
---|
554 | ! WHERE(tot_7doc(ig,lg,:,:).GT.zero) |
---|
555 | ! DOC(ig,lg,jg,:,:,icarbon)=doc_old(ig,lg,jg,:,:)-erodoc(ig,lg,:,:)*doc_old(ig,lg,jg,:,:)/tot_7doc(ig,lg,:,:) + & |
---|
556 | ! & (erodep(ig,lg)/dly(8)* doc_old(ig,lg,8,:,:))*doc_old(ig,lg,jg,:,:)/tot_7doc(ig,lg,:,:) |
---|
557 | ! ENDWHERE |
---|
558 | |
---|
559 | ! WRITE(numout,*) 'SUMero3',ig,lg,SUM(litter_below(ig,:,lg,:,icarbon)),SUM(carbon(ig,:,lg,jg)), & |
---|
560 | ! lignin_struc_below(ig,lg,jg),SUM(DOC(ig,lg,jg,:,:,icarbon)) |
---|
561 | ENDDO ! top 7 soil layers |
---|
562 | DO jg=nslm-1,8,-1 ! top 8-10 soil layers |
---|
563 | carbon(ig,:,lg,jg)= carbon_old(ig,:,lg,jg)-erodep(ig,lg)/dly(jg)*carbon_old(ig,:,lg,jg)+ & |
---|
564 | & erodep(ig,lg)/dly(jg+1)*carbon_old(ig,:,lg,jg+1) |
---|
565 | litter_below(ig,:,lg,jg,icarbon)=litbe_old(ig,:,lg,jg)-erodep(ig,lg)/dly(jg)*litbe_old(ig,:,lg,jg) + & |
---|
566 | & erodep(ig,lg)/dly(jg+1)*litbe_old(ig,:,lg,jg+1) |
---|
567 | ! IF (SUM(litter_below(ig,:,lg,jg,icarbon)).GT.zero) THEN |
---|
568 | ! lignin_struc_below(ig,lg,jg)= (lignin_struc_below(ig,lg,jg)*(litbe_tot(ig,lg,jg)-erodep(ig,lg)/dly(jg)*litbe_tot(ig,lg,jg)) + & |
---|
569 | ! & lignin_struc_below(ig,lg,jg+1)*erodep(ig,lg)/dly(jg+1)*litbe_tot(ig,lg,jg+1)) / SUM(litter_below(ig,:,lg,jg,icarbon)) |
---|
570 | ! ENDIF |
---|
571 | ! DOC(ig,lg,jg,:,:,icarbon)= doc_old(ig,lg,jg,:,:)-erodep(ig,lg)/dly(jg)*doc_old(ig,lg,jg,:,:) + & |
---|
572 | ! & erodep(ig,lg)/dly(jg+1)*doc_old(ig,lg,jg+1,:,:) |
---|
573 | |
---|
574 | ! WRITE(numout,*) 'SUMero4',ig,lg,SUM(litter_below(ig,:,lg,:,icarbon)),SUM(carbon(ig,:,lg,jg)), & |
---|
575 | ! lignin_struc_below(ig,lg,jg),SUM(DOC(ig,lg,jg,:,:,icarbon)) |
---|
576 | ENDDO |
---|
577 | ! |
---|
578 | IF (depth_deepsoil(ig,lg).GT.min_sechiba .AND.depth_deepsoil(ig,lg).GT.erodep(ig,lg)) THEN |
---|
579 | carbon(ig,:,lg,nslm)=carbon_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*carbon_old(ig,:,lg,nslm) + & |
---|
580 | & erodep(ig,lg)/depth_deepsoil(ig,lg)*carbon_old(ig,:,lg,nslmd) |
---|
581 | carbon(ig,:,lg,nslmd)=(1.0-erodep(ig,lg)/depth_deepsoil(ig,lg))*carbon_old(ig,:,lg,nslmd) |
---|
582 | |
---|
583 | litter_below(ig,:,lg,nslm,icarbon)=litbe_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*litbe_old(ig,:,lg,nslm) + & |
---|
584 | & erodep(ig,lg)/depth_deepsoil(ig,lg)*litbe_old(ig,:,lg,nslmd) |
---|
585 | litter_below(ig,:,lg,nslmd,icarbon)=(1.0-erodep(ig,lg)/depth_deepsoil(ig,lg))*litbe_old(ig,:,lg,nslmd) |
---|
586 | IF (SUM(litter_below(ig,:,lg,nslm,icarbon)).GT.zero) THEN |
---|
587 | ! lignin_struc_below(ig,lg,nslm)=(lignin_struc_below(ig,lg,nslm)*(litbe_tot(ig,lg,nslm)-erodep(ig,lg)/dly(nslm)* & |
---|
588 | ! & litbe_tot(ig,lg,nslm))+erodep(ig,lg)/depth_deepsoil(ig,lg)*litbe_tot(ig,lg,nslmd)* & |
---|
589 | ! & lignin_struc_below(ig,lg,nslmd))/SUM(litter_below(ig,:,lg,nslm,icarbon)) |
---|
590 | ENDIF |
---|
591 | ! DOC(ig,lg,nslm,:,:,icarbon)=doc_old(ig,lg,nslm,:,:)-erodep(ig,lg)/dly(nslm)*doc_old(ig,lg,nslm,:,:) + & |
---|
592 | ! & erodep(ig,lg)/depth_deepsoil(ig,lg)*doc_old(ig,lg,nslmd,:,:) |
---|
593 | ! DOC(ig,lg,nslmd,:,:,icarbon)=(1.0-erodep(ig,lg)/depth_deepsoil(ig,lg))*doc_old(ig,lg,nslmd,:,:) |
---|
594 | |
---|
595 | ! WRITE(numout,*) 'SUMero5',ig,lg,SUM(litter_below(ig,:,lg,nslm:nslmd,icarbon)),SUM(carbon(ig,:,lg,nslm:nslmd)), & |
---|
596 | ! lignin_struc_below(ig,lg,nslm),lignin_struc_below(ig,lg,nslmd),SUM(DOC(ig,lg,nslm:nslmd,:,:,icarbon)) |
---|
597 | ! kjpindex,nvm,nslmd,ndoc,npool,nelements |
---|
598 | ELSEIF (depth_deepsoil(ig,lg).LE.erodep(ig,lg)) THEN |
---|
599 | carbon(ig,:,lg,nslm)=carbon_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*carbon_old(ig,:,lg,nslm) + & |
---|
600 | & carbon_old(ig,:,lg,nslmd) |
---|
601 | carbon(ig,:,lg,nslmd)=zero |
---|
602 | |
---|
603 | litter_below(ig,:,lg,nslm,icarbon)=litbe_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*litbe_old(ig,:,lg,nslm) + & |
---|
604 | & litbe_old(ig,:,lg,nslmd) |
---|
605 | litter_below(ig,:,lg,nslmd,icarbon)= zero |
---|
606 | ! IF (SUM(litter_below(ig,:,lg,jg,icarbon)).GT.zero) THEN |
---|
607 | ! lignin_struc_below(ig,lg,nslm)=(lignin_struc_below(ig,lg,nslm)*(litbe_tot(ig,lg,nslm)-erodep(ig,lg)/dly(nslm)* & |
---|
608 | ! & litbe_tot(ig,lg,nslm))+litbe_tot(ig,lg,nslmd)*lignin_struc_below(ig,lg,nslmd))/ & |
---|
609 | ! & SUM(litter_below(ig,:,lg,nslm,icarbon)) |
---|
610 | ! ENDIF |
---|
611 | lignin_struc_below(ig,lg,nslmd)=zero |
---|
612 | |
---|
613 | ! DOC(ig,lg,nslm,:,:,icarbon)=doc_old(ig,lg,nslm,:,:)-erodep(ig,lg)/dly(nslm)*doc_old(ig,lg,nslm,:,:) + & |
---|
614 | ! & doc_old(ig,lg,nslmd,:,:) |
---|
615 | ! DOC(ig,lg,nslmd,:,:,icarbon)=zero |
---|
616 | |
---|
617 | ! WRITE(numout,*) 'SUMero6',ig,lg,SUM(litter_below(ig,:,lg,nslm:nslmd,icarbon)),SUM(carbon(ig,:,lg,nslm:nslmd)), & |
---|
618 | ! lignin_struc_below(ig,lg,nslm),lignin_struc_below(ig,lg,nslmd),SUM(DOC(ig,lg,nslm:nslmd,:,:,icarbon)) |
---|
619 | ELSE |
---|
620 | carbon(ig,:,lg,nslm)=carbon_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*carbon_old(ig,:,lg,nslm) |
---|
621 | carbon(ig,:,lg,nslmd)=zero |
---|
622 | |
---|
623 | litter_below(ig,:,lg,nslm,icarbon)=litbe_old(ig,:,lg,nslm)-erodep(ig,lg)/dly(nslm)*litbe_old(ig,:,lg,nslm) |
---|
624 | litter_below(ig,:,lg,nslmd,icarbon)= zero |
---|
625 | lignin_struc_below(ig,lg,nslmd)=zero |
---|
626 | |
---|
627 | ! DOC(ig,lg,nslm,:,:,icarbon)=doc_old(ig,lg,nslm,:,:)-erodep(ig,lg)/dly(nslm)*doc_old(ig,lg,nslm,:,:) |
---|
628 | ! DOC(ig,lg,nslmd,:,:,icarbon)=zero |
---|
629 | |
---|
630 | ! WRITE(numout,*) 'SUMero7',ig,lg,SUM(litter_below(ig,:,lg,nslm:nslmd,icarbon)),SUM(carbon(ig,:,lg,nslm:nslmd)), & |
---|
631 | ! lignin_struc_below(ig,lg,nslm),lignin_struc_below(ig,lg,nslmd),SUM(DOC(ig,lg,nslm:nslmd,:,:,icarbon)) |
---|
632 | ENDIF |
---|
633 | depth_deepsoil(ig,lg)=MAX(zero,depth_deepsoil(ig,lg)-erodep(ig,lg)) |
---|
634 | |
---|
635 | ! WRITE(numout,*) 'soilcarbon_passive', carbon(ig,3,lg,1:7) |
---|
636 | ! WRITE(numout,*) 'soilcarbon_active', carbon(ig,1,lg,1:7) |
---|
637 | ENDIF |
---|
638 | !WRITE(numout,*) 'ero_totc', ig,lg,SUM(litter_below(ig,:,lg,:,icarbon)),SUM(litter_below(ig,:,1,:,icarbon)), & |
---|
639 | !& SUM(lignin_struc_below(ig,lg,:)) |
---|
640 | ! WRITE(numout,*) 'ero_totc', veget_max(ig,lg),MinVal(carbon(ig,:,lg,1:7)),SUM(tot_eroc(ig,:)) |
---|
641 | ENDDO !DO lg=1,nvm |
---|
642 | ! WRITE(numout,*) 'erosion carbon', SUM(tot_eroc(ig,:)) |
---|
643 | ENDDO !DO ig=1,nbpt |
---|
644 | eroc_active(:,:)=eroc(:,iactive,:) |
---|
645 | eroc_slow(:,:)=eroc(:,islow,:) |
---|
646 | eroc_passive(:,:)=eroc(:,ipassive,:) |
---|
647 | |
---|
648 | !WRITE(numout,*) 'SUMero2_litter',MAXVAL(litter_above),MAXVAL(litter_below) |
---|
649 | !WRITE(numout,*) 'SUMero2_SOC',MAXVAL(carbon) |
---|
650 | !WRITE(numout,*) 'SUMero2_DOC',MAXVAL(DOC) |
---|
651 | |
---|
652 | ! 2.5 Average export rate of sediment and POC from land to river network |
---|
653 | SED_EXP_agg(:,:)=zero |
---|
654 | POC_EXP_agg(:,:)=zero |
---|
655 | DOC_ERO_agg(:,:)=zero |
---|
656 | sed_text(:,:)=zero |
---|
657 | |
---|
658 | ! *** Assumption-1: the enrichment ratio of clay fraction is enrich_clay=3 |
---|
659 | ! In case the sum of fractions of clay, silt, sand are less than un |
---|
660 | sed_text(:,ic_clay)=MAX(zero,un-textfrac(:,ic_silt)-textfrac(:,ic_sand)) |
---|
661 | |
---|
662 | sed_text(:,ic_clay)=MAX(sed_text(:,ic_clay),MIN(un,sed_text(:,ic_clay)*enrich_clay)) |
---|
663 | sed_text(:,ic_silt)=MIN(textfrac(:,ic_silt)*enrich_silt,un-sed_text(:,ic_clay)) |
---|
664 | sed_text(:,ic_sand)=MIN(textfrac(:,ic_sand),un-sed_text(:,ic_clay)-sed_text(:,ic_silt)) |
---|
665 | |
---|
666 | ! *** Assumption-2: All sediment transported into river channel is the clay fraction; Or we just simulate the |
---|
667 | ! total sediment transport, without distinguish the clay, silt and sand fraction |
---|
668 | ! sed_text(:,ic_clay)= un |
---|
669 | ! sed_text(:,ic_silt)= zero |
---|
670 | ! sed_text(:,ic_sand)= zero |
---|
671 | |
---|
672 | DO lg=1,nvm |
---|
673 | !WHERE(erorate(:,lg).GT.zero) |
---|
674 | SED_EXP_agg(:,ic_clay)=SED_EXP_agg(:,ic_clay)+ & |
---|
675 | erorate(:,lg)*sed_text(:,ic_clay)*veget_max(:,lg) |
---|
676 | |
---|
677 | SED_EXP_agg(:,ic_silt)=SED_EXP_agg(:,ic_silt)+ & |
---|
678 | erorate(:,lg)*sed_text(:,ic_silt)*veget_max(:,lg) |
---|
679 | |
---|
680 | SED_EXP_agg(:,ic_sand)=SED_EXP_agg(:,ic_sand)+ & |
---|
681 | erorate(:,lg)*sed_text(:,ic_sand)*veget_max(:,lg) |
---|
682 | !ENDWHERE |
---|
683 | ! |
---|
684 | DO icarb=1,ncarb |
---|
685 | ! WHERE(eroc(:,icarb,lg).GT.zero) |
---|
686 | POC_EXP_agg(:,icarb)=POC_EXP_agg(:,icarb)+eroc(:,icarb,lg)*veget_max(:,lg) |
---|
687 | ! ENDWHERE |
---|
688 | ENDDO |
---|
689 | ! |
---|
690 | DO ii = 1,ndoc |
---|
691 | DOC_ERO_agg(:,iactive) = DOC_ERO_agg(:,iactive)+ & |
---|
692 | (erodoc(:,lg,ii,imetabo)+erodoc(:,lg,ii,imetbel)+erodoc(:,lg,ii,iact))*veget_max(:,lg) |
---|
693 | DOC_ERO_agg(:,islow) = DOC_ERO_agg(:,islow)+ & |
---|
694 | (erodoc(:,lg,ii,istrabo)+erodoc(:,lg,ii,istrbel)+erodoc(:,lg,ii,islo))*veget_max(:,lg) |
---|
695 | DOC_ERO_agg(:,ipassive) = DOC_ERO_agg(:,ipassive)+ & |
---|
696 | erodoc(:,lg,ii,ipas)*veget_max(:,lg) |
---|
697 | ENDDO |
---|
698 | ! |
---|
699 | ENDDO !DO lg=1,nvm |
---|
700 | !SED_EXP_agg(:,:)=SED_EXP_agg(:,:) |
---|
701 | !POC_EXP_agg(:,:)=POC_EXP_agg(:,:) |
---|
702 | !DOC_ERO_agg(:,:)=DOC_ERO_agg(:,:) |
---|
703 | |
---|
704 | !WRITE(numout,*) 'ZHCero_Cexp0' |
---|
705 | !DO ig=1,nbpt |
---|
706 | ! WRITE(numout,*) 'Runoff', ig,totrunoff_d(ig) |
---|
707 | ! WRITE(numout,*) 'effgrid_perc',effgrid_perc(ig) |
---|
708 | ! WRITE(numout,*) 'std_totsed',std_totsed(ig) |
---|
709 | ! WRITE(numout,*) 'erocap_pft',ig,SUM(erocap(ig,:)) |
---|
710 | ! WRITE(numout,*) 'effgrid_perc',effgrid_perc(ig) |
---|
711 | ! WRITE(numout,*) 'Erorate', erorate(ig,:) |
---|
712 | ! WRITE(numout,*) 'textfrac',textfrac(ig,:),sed_text(ig,:) |
---|
713 | ! WRITE(numout,*) 'Sediment', SED_EXP_agg(ig,:) |
---|
714 | ! WRITE(numout,*) 'POC', POC_EXP_agg(ig,:)*one_day/dt_sechiba |
---|
715 | ! WRITE(numout,*) 'Concent_g/m^3', SED_EXP_agg(ig,:)*one_day/dt_sechiba/totrunoff_d(ig)*1000.0 |
---|
716 | !ENDDO |
---|
717 | !! 2.6 Update SOC, litter_below and DOC pools caused by sediment deposition |
---|
718 | seddep_rate(:,:)=zero |
---|
719 | soc_deposition(:,:,:)=zero |
---|
720 | |
---|
721 | sed_deposition_tot(:) = (sed_deposition_d(:,ic_clay)+sed_deposition_d(:,ic_silt)+ & |
---|
722 | & sed_deposition_d(:,ic_sand))*dt_erosion/one_day |
---|
723 | poc_deposition_d(:,:) = poc_deposition_d(:,:)*dt_erosion/one_day |
---|
724 | |
---|
725 | max_depth=min(laydep(7),dly(8)) |
---|
726 | DO lg=2,nvm |
---|
727 | WHERE (bio_frac(:).GT.zero.AND.bulkdens(:).GT.zero) |
---|
728 | seddep_rate(:,lg) = sed_deposition_tot(:)/bio_frac(:)/bulkdens(:) |
---|
729 | ENDWHERE |
---|
730 | ! H.Zhang: We assumed that 1) the deposited sediemnt and POC will evenly distributed in |
---|
731 | ! each grid-cell, 2) the sediment and POC are fully mixed; 3) the depth of daily sediment |
---|
732 | ! deopisiton does not exceed laydep(4) (~5.0 cm) |
---|
733 | ! This assumption did affect the C balance |
---|
734 | WHERE (seddep_rate(:,lg).GT.max_depth) |
---|
735 | seddep_rate(:,lg) = max_depth-min_sechiba |
---|
736 | ENDWHERE |
---|
737 | depth_deepsoil(:,lg)=depth_deepsoil(:,lg)+seddep_rate(:,lg) ! unit: m |
---|
738 | ! |
---|
739 | DO icarb=1,ncarb |
---|
740 | WHERE (bio_frac(:).GT.zero.AND.veget_max(:,lg).GT.zero) |
---|
741 | soc_deposition(:,icarb,lg)=1000.*poc_deposition_d(:,icarb)/bio_frac(:) ! Unit: g C m-2 day-1 |
---|
742 | ENDWHERE |
---|
743 | ENDDO !icarb=1,ncarb |
---|
744 | ENDDO !! DO lg=2,nvm |
---|
745 | |
---|
746 | carbon_old(:,:,:,:)=carbon(:,:,:,:) |
---|
747 | litbe_old(:,:,:,:)=litter_below(:,:,:,:,icarbon) |
---|
748 | litbe_tot(:,:,:)=SUM(litbe_old(:,:,:,:),DIM=2) |
---|
749 | doc_old(:,:,:,:,:)=doc(:,:,:,:,:,icarbon) |
---|
750 | tot_7c(:,:,:)=zero |
---|
751 | tot_7litbelow(:,:,:)=zero |
---|
752 | tot_7doc(:,:,:,:)=zero |
---|
753 | |
---|
754 | DO ig=1,nbpt |
---|
755 | DO lg=2,nvm |
---|
756 | ! |
---|
757 | !WRITE(numout,*) 'ZHC_below0', ig,lg,seddep_rate(ig,lg) |
---|
758 | IF ((veget_max(ig,lg) .GT. 0).AND.(seddep_rate(ig,lg).GT.zero).AND.(MinVal(SUM(carbon(ig,:,lg,1:7),DIM=2)).GT.0)) THEN !!! 1<=fldir<=128 |
---|
759 | DO jg=1,7 |
---|
760 | tot_7c(ig,:,lg)=tot_7c(ig,:,lg)+carbon(ig,:,lg,jg) |
---|
761 | tot_7litbelow(ig,:,lg)=tot_7litbelow(ig,:,lg)+litter_below(ig,:,lg,jg,icarbon) |
---|
762 | tot_7doc(ig,lg,:,:)=tot_7doc(ig,lg,:,:)+DOC(ig,lg,jg,:,:,icarbon) |
---|
763 | ENDDO !DO jg=,7 |
---|
764 | ! Situation-1: Part of original layer-1 enter layer-2; newly deposited soil does not reach layer-2 |
---|
765 | carbon(ig,:,lg,nslmd)=carbon_old(ig,:,lg,nslmd)+seddep_rate(ig,lg)/dly(nslm)*carbon_old(ig,:,lg,nslm) |
---|
766 | litter_below(ig,:,lg,nslmd,icarbon)=litbe_old(ig,:,lg,nslmd)+seddep_rate(ig,lg)/dly(nslm)*litbe_old(ig,:,lg,nslm) |
---|
767 | IF (SUM(litter_below(ig,:,lg,nslmd,icarbon)).GT.zero) THEN |
---|
768 | lignin_struc_below(ig,lg,nslmd)=(litbe_tot(ig,lg,nslmd)*lignin_struc_below(ig,lg,nslmd)+ & |
---|
769 | & seddep_rate(ig,lg)/dly(nslm)*litbe_tot(ig,lg,nslm)*lignin_struc_below(ig,lg,nslm))/ & |
---|
770 | & SUM(litter_below(ig,:,lg,nslmd,icarbon)) |
---|
771 | ENDIF |
---|
772 | DOC(ig,lg,nslmd,:,:,icarbon)=doc_old(ig,lg,nslmd,:,:)+seddep_rate(ig,lg)/dly(nslm)*doc_old(ig,lg,nslm,:,:) |
---|
773 | ! |
---|
774 | DO jg=nslm,9,-1 |
---|
775 | carbon(ig,:,lg,jg)=(dly(jg)-seddep_rate(ig,lg))/dly(jg)*carbon_old(ig,:,lg,jg)+seddep_rate(ig,lg)/laydep(jg-1)*carbon_old(ig,:,lg,jg-1) |
---|
776 | litter_below(ig,:,lg,jg,icarbon)=(dly(jg)-seddep_rate(ig,lg))/dly(jg)*litbe_old(ig,:,lg,jg)+seddep_rate(ig,lg)/laydep(jg-1)*litbe_old(ig,:,lg,jg-1) |
---|
777 | ! IF (SUM(litter_below(ig,:,lg,jg,icarbon)).GT.zero) THEN |
---|
778 | ! lignin_struc_below(ig,lg,jg)=(lignin_struc_below(ig,lg,jg)*(dly(jg)-seddep_rate(ig,lg))/dly(jg)*litbe_tot(ig,lg,jg)+ & |
---|
779 | ! & lignin_struc_below(ig,lg,jg-1)*seddep_rate(ig,lg)/laydep(jg-1)*litbe_tot(ig,lg,jg-1))/SUM(litter_below(ig,:,lg,jg,icarbon)) |
---|
780 | ! ENDIF |
---|
781 | DOC(ig,lg,jg,:,:,icarbon)=(dly(jg)-seddep_rate(ig,lg))/dly(jg)*doc_old(ig,lg,jg,:,:)+seddep_rate(ig,lg)/laydep(jg-1)*doc_old(ig,lg,jg-1,:,:) |
---|
782 | ENDDO ! DO jg=nslm,2,-1 |
---|
783 | carbon(ig,:,lg,8)=(dly(8)-seddep_rate(ig,lg))/dly(8)*carbon_old(ig,:,lg,8)+seddep_rate(ig,lg)/laydep(7)*tot_7c(ig,:,lg) |
---|
784 | litter_below(ig,:,lg,8,icarbon)=(dly(8)-seddep_rate(ig,lg))/dly(8)*litbe_old(ig,:,lg,8)+seddep_rate(ig,lg)/laydep(7)*tot_7litbelow(ig,:,lg) |
---|
785 | DOC(ig,lg,8,:,:,icarbon)=(dly(8)-seddep_rate(ig,lg))/dly(8)*doc_old(ig,lg,8,:,:)+seddep_rate(ig,lg)/laydep(7)*tot_7doc(ig,lg,:,:) |
---|
786 | ! |
---|
787 | DO jg=7,1,-1 |
---|
788 | WHERE (tot_7c(ig,:,lg).GT.zero) |
---|
789 | carbon(ig,:,lg,1) = ((un-seddep_rate(ig,lg)/laydep(7))*tot_7c(ig,:,lg)+soc_deposition(ig,:,lg))*(carbon_old(ig,:,lg,jg)/tot_7c(ig,:,lg)) |
---|
790 | ENDWHERE |
---|
791 | litter_below(ig,:,lg,jg,icarbon) = (un-seddep_rate(ig,lg)/laydep(7))*litbe_old(ig,:,lg,jg) ! No litterfall depositioon |
---|
792 | DOC(ig,lg,jg,:,:,icarbon) = (un-seddep_rate(ig,lg)/laydep(7))*doc_old(ig,lg,jg,:,:) ! DOC will enter into the nslm soil layer following water flow |
---|
793 | ENDDO |
---|
794 | !WRITE(numout,*) 'ZHC_below10', ig,lg,SUM(litter_below(ig,:,lg,nslm,icarbon)),SUM(lignin_struc_below(ig,lg,:)) |
---|
795 | ENDIF |
---|
796 | |
---|
797 | ! Make sure the lignin_struc_ is smaller than 1.0 and higher than 0.0 |
---|
798 | lignin_struc_above(ig,lg)=MIN(MAX(lignin_struc_above(ig,lg),zero),1.0) |
---|
799 | DO jg=nslmd,1,-1 |
---|
800 | lignin_struc_below(ig,lg,jg) = MIN(MAX(lignin_struc_below(ig,lg,jg),zero),1.0) |
---|
801 | ENDDO |
---|
802 | |
---|
803 | ENDDO !! DO lg=2,nvm |
---|
804 | ENDDO !DO ig=1,nbpt |
---|
805 | |
---|
806 | !WRITE(numout,*) 'SUMero3_litter',SUM(litter_above),SUM(litter_below),SUM(lignin_struc_below) |
---|
807 | !WRITE(numout,*) 'SUMero3_SOC',SUM(carbon) |
---|
808 | !WRITE(numout,*) 'SUMero3_DOC',SUM(DOC) |
---|
809 | !! 2.7 Write soil & carbon erosion rate to the history file |
---|
810 | |
---|
811 | CALL xios_orchidee_send_field("bulkdens",bulkdens*one_day/dt_erosion) ! dt_erosion=one_day |
---|
812 | CALL xios_orchidee_send_field("totrunoff_d",totrunoff_d*one_day/dt_erosion) |
---|
813 | CALL xios_orchidee_send_field("peakrunoff_d",peakrunoff_d) |
---|
814 | CALL xios_orchidee_send_field("totrundra_d",totrundra_d*one_day/dt_erosion) |
---|
815 | CALL xios_orchidee_send_field("peakrundra_d",peakrundra_d) |
---|
816 | CALL xios_orchidee_send_field("cf_veget",cf_veget) |
---|
817 | CALL xios_orchidee_send_field("kf_litter",kf_litter) |
---|
818 | CALL xios_orchidee_send_field("kf_root",kf_root) |
---|
819 | CALL xios_orchidee_send_field("erocap",erocap*one_day/dt_erosion) |
---|
820 | CALL xios_orchidee_send_field("erodep",erodep*one_day/dt_erosion) |
---|
821 | CALL xios_orchidee_send_field("eroc_active",eroc_active*one_day/dt_erosion) |
---|
822 | CALL xios_orchidee_send_field("eroc_slow",eroc_slow*one_day/dt_erosion) |
---|
823 | CALL xios_orchidee_send_field("eroc_passive",eroc_passive*one_day/dt_erosion) |
---|
824 | CALL xios_orchidee_send_field("tot_erocap",tot_erocap*one_day/dt_erosion) |
---|
825 | CALL xios_orchidee_send_field("ave_erodep",ave_erodep*one_day/dt_erosion) |
---|
826 | CALL xios_orchidee_send_field("tot_eroc",tot_eroc*one_day/dt_erosion) |
---|
827 | CALL xios_orchidee_send_field("depth_deepsoil",depth_deepsoil*one_day/dt_erosion) |
---|
828 | CALL xios_orchidee_send_field("seddep_rate",seddep_rate*one_day/dt_erosion) |
---|
829 | ! |
---|
830 | ! Reinitialize values |
---|
831 | time_counter_ero = zero |
---|
832 | cveg_f = zero |
---|
833 | cf_veget=zero |
---|
834 | kflitter=zero |
---|
835 | kfroot=zero |
---|
836 | kf_litter = zero |
---|
837 | kf_root = zero |
---|
838 | peakrunoff = zero |
---|
839 | meanrunoff = zero |
---|
840 | peakrundra = zero |
---|
841 | meanrundra = zero |
---|
842 | erocap = zero |
---|
843 | erodep = zero |
---|
844 | eroc=zero |
---|
845 | tot_erocap=zero |
---|
846 | ave_erodep =zero |
---|
847 | tot_eroc =zero |
---|
848 | sed_deposition_d = zero |
---|
849 | sed_deposition_tot= zero |
---|
850 | poc_deposition_d = zero |
---|
851 | no_events = 0 |
---|
852 | |
---|
853 | ! WRITE(numout,*) 'Erosion_ZHC100' |
---|
854 | ! WRITE(numout,*) 'litt_ab: ', MAXVAL(litter_above),MINVAL(litter_above) |
---|
855 | ! WRITE(numout,*) 'litt_be: ', MAXVAL(litter_below),MINVAL(litter_below) |
---|
856 | ! WRITE(numout,*) 'ligninbe',MAXVAL(lignin_struc_below),MINVAL(lignin_struc_below) |
---|
857 | ! WRITE(numout,*) 'DOC',MAXVAL(DOC),MINVAL(DOC) |
---|
858 | ! WRITE(numout,*) 'SOC',MAXVAL(carbon),MINVAL(carbon) |
---|
859 | ! |
---|
860 | ! WRITE(numout,*) 'End of this day' |
---|
861 | ENDIF !IF NINT(time_counter_ero) .GE. NINT(dt_routing) |
---|
862 | |
---|
863 | ! |
---|
864 | END SUBROUTINE erosion_main |
---|
865 | ! |
---|
866 | ! |
---|
867 | !! ================================================================================================================================ |
---|
868 | !! SUBROUTINE : erosion_clear |
---|
869 | !! |
---|
870 | !> BRIEF : This subroutine deallocates the block memory previously allocated. |
---|
871 | !! \n |
---|
872 | !_ ================================================================================================================================ |
---|
873 | SUBROUTINE erosion_clear() |
---|
874 | l_first_erosion = .TRUE. |
---|
875 | ! |
---|
876 | IF (ALLOCATED(gridarea)) DEALLOCATE(gridarea) |
---|
877 | IF (ALLOCATED(cveg_f)) DEALLOCATE(cveg_f) |
---|
878 | IF (ALLOCATED(kflitter)) DEALLOCATE(kflitter) |
---|
879 | IF (ALLOCATED(kfroot)) DEALLOCATE(kfroot) |
---|
880 | IF (ALLOCATED(peakrunoff)) DEALLOCATE(peakrunoff) |
---|
881 | IF (ALLOCATED(meanrunoff)) DEALLOCATE(meanrunoff) |
---|
882 | IF (ALLOCATED(peakrundra)) DEALLOCATE(peakrundra) |
---|
883 | IF (ALLOCATED(meanrundra)) DEALLOCATE(meanrundra) |
---|
884 | IF (ALLOCATED(no_events)) DEALLOCATE(no_events) |
---|
885 | IF (ALLOCATED(eroc_active)) DEALLOCATE(eroc_active) |
---|
886 | IF (ALLOCATED(eroc_slow)) DEALLOCATE(eroc_slow) |
---|
887 | IF (ALLOCATED(eroc_passive)) DEALLOCATE(eroc_passive) |
---|
888 | IF (ALLOCATED(tot_7c)) DEALLOCATE(tot_7c) |
---|
889 | IF (ALLOCATED(tot_7litbelow)) DEALLOCATE(tot_7litbelow) |
---|
890 | IF (ALLOCATED(tot_7doc)) DEALLOCATE(tot_7doc) |
---|
891 | !**************** output variable ************************* |
---|
892 | END SUBROUTINE erosion_clear |
---|
893 | ! |
---|
894 | ! |
---|
895 | ! |
---|
896 | !! ================================================================================================================================ |
---|
897 | !! SUBROUTINE : erosionmap |
---|
898 | !! |
---|
899 | !> BRIEF : This subroutine reads and interpolates erosion map |
---|
900 | !! \n |
---|
901 | !_ ================================================================================================================================ |
---|
902 | SUBROUTINE erosionmap(nbpt,index,lalo,neighbours,resolution,contfrac,std_totsed,gridarea,effgrid_perc) |
---|
903 | ! |
---|
904 | IMPLICIT NONE |
---|
905 | ! |
---|
906 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
907 | INTEGER(i_std), INTENT(in) :: index(nbpt) !! Index on the global map. |
---|
908 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) |
---|
909 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
910 | ! INTEGER(i_std), INTENT(out) :: maxnbp !! maxnbp(1):number of fine grid in a ORCHIDEE grid, |
---|
911 | !! 2,3: total row and column number of sub_grid |
---|
912 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) |
---|
913 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
914 | CHARACTER(LEN=80) :: filename !! Name of the netcdf file (unitless) |
---|
915 | ! INTEGER(i_std) :: nbpmax, nix, njx, fopt !! IDices, (unitless) |
---|
916 | CHARACTER(LEN=30) :: callsign !! |
---|
917 | REAL(r_std), INTENT(out) :: std_totsed(nbpt) !! standard sediment discharge amount rate for each ORCHIDEE cell(ton) |
---|
918 | REAL(r_std), INTENT(out) :: effgrid_perc(nbpt) !! percentage of the effective area in a ORCHIDEE pixel (0-1) |
---|
919 | REAL(r_std), INTENT(out) :: gridarea(nbpt) !! area of each grid cell (m^2) |
---|
920 | ! |
---|
921 | INTEGER :: ALLOC_ERR !! |
---|
922 | LOGICAL :: ok_interpol = .FALSE. !! Flag for interpolation (true/false) |
---|
923 | LOGICAL :: is_root_prc_old !! |
---|
924 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib,ig, jb, ip, jp, itype,subres !! Indices (unitless) |
---|
925 | REAL(r_std) :: lev(1), date, dt, coslat !! |
---|
926 | REAL(r_std) :: lon0,lat0,res_lon,res_lat,cellarea !! size (m) of target grid at longitude and latitude direction |
---|
927 | REAL(r_std) :: sumsed !! Total erosion from one ORCHIDEE cell in the simulation resolution (ton) |
---|
928 | REAL(r_std) :: sumarea !! total effective area (m^2) |
---|
929 | INTEGER(i_std) :: itau(1) !! |
---|
930 | |
---|
931 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: latrel !! Latitude |
---|
932 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lonrel !! Longitude |
---|
933 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu !! Resolution read on the map |
---|
934 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: totsed !! erosion rate for the finer grid, read from forcing data(ton) |
---|
935 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: effgrid !! percentage of the effective area in a ORCHIDEE pixel (0-1) |
---|
936 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask !! Mask to exclude some points (unitless) |
---|
937 | !sub_gridarea, sub_index |
---|
938 | !_ ================================================================================================================================ |
---|
939 | ! |
---|
940 | is_root_prc_old=is_root_prc |
---|
941 | is_root_prc=.TRUE. |
---|
942 | ! |
---|
943 | filename = 'erosion.nc' ! |
---|
944 | CALL getin_p('EROSION_FILE',filename) |
---|
945 | ! |
---|
946 | WRITE(numout,*) 'read info of erosion map', nbpt,filename |
---|
947 | IF (is_root_prc) THEN |
---|
948 | CALL flininfo(filename,iml, jml, lml, tml, fid) |
---|
949 | CALL flinclo(fid) |
---|
950 | ENDIF |
---|
951 | CALL bcast(iml) |
---|
952 | CALL bcast(jml) |
---|
953 | CALL bcast(lml) |
---|
954 | CALL bcast(tml) |
---|
955 | WRITE(numout,*) 'read info of erosion map', nbpt,is_root_prc, filename |
---|
956 | ! |
---|
957 | ALLOC_ERR=-1 |
---|
958 | ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR) |
---|
959 | IF (ALLOC_ERR/=0) THEN |
---|
960 | WRITE(numout,*) "ERROR IN ALLOCATION of lonrel: ",ALLOC_ERR, jml |
---|
961 | STOP |
---|
962 | ENDIF |
---|
963 | lonrel(:,:)=zero |
---|
964 | ! |
---|
965 | ALLOC_ERR=-1 |
---|
966 | ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR) |
---|
967 | IF (ALLOC_ERR/=0) THEN |
---|
968 | WRITE(numout,*) "ERROR IN ALLOCATION of latrel: ",ALLOC_ERR, jml |
---|
969 | STOP |
---|
970 | ENDIF |
---|
971 | latrel(:,:)=zero |
---|
972 | ! |
---|
973 | ALLOC_ERR=-1 |
---|
974 | ALLOCATE (totsed(iml,jml), STAT=ALLOC_ERR) |
---|
975 | IF (ALLOC_ERR/=0) THEN |
---|
976 | WRITE(numout,*) "ERROR IN ALLOCATION of totsed: ",ALLOC_ERR, jml |
---|
977 | STOP |
---|
978 | ENDIF |
---|
979 | totsed(:,:)=zero |
---|
980 | ! |
---|
981 | ALLOC_ERR=-1 |
---|
982 | ALLOCATE (effgrid(iml,jml), STAT=ALLOC_ERR) |
---|
983 | IF (ALLOC_ERR/=0) THEN |
---|
984 | WRITE(numout,*) "ERROR IN ALLOCATION of effgrid: ",ALLOC_ERR, jml |
---|
985 | STOP |
---|
986 | ENDIF |
---|
987 | effgrid(:,:)=zero |
---|
988 | ! |
---|
989 | ALLOC_ERR=-1 |
---|
990 | ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR) |
---|
991 | IF (ALLOC_ERR/=0) THEN |
---|
992 | WRITE(numout,*) "ERROR IN ALLOCATION of resol_lu : ",ALLOC_ERR |
---|
993 | STOP |
---|
994 | ENDIF |
---|
995 | resol_lu(:,:,:)=zero |
---|
996 | ! |
---|
997 | ALLOC_ERR=-1 |
---|
998 | ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) |
---|
999 | IF (ALLOC_ERR/=0) THEN |
---|
1000 | WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR |
---|
1001 | STOP |
---|
1002 | ENDIF |
---|
1003 | mask(:,:) = 0. |
---|
1004 | ! |
---|
1005 | ! ! |
---|
1006 | WRITE(numout,*) "Open erosion file",is_root_prc, fid |
---|
1007 | IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid) |
---|
1008 | WRITE(numout,*) "Open erosion file1",iml, jml, lml,fid |
---|
1009 | CALL bcast(lonrel) |
---|
1010 | CALL bcast(latrel) |
---|
1011 | CALL bcast(lev) |
---|
1012 | CALL bcast(itau) |
---|
1013 | CALL bcast(date) |
---|
1014 | CALL bcast(dt) |
---|
1015 | ! |
---|
1016 | !WRITE(numout,*) "Open erosion totsed",iml, jml, lml, tml,is_root_prc,fid |
---|
1017 | IF (is_root_prc) CALL flinget(fid, 'totsed', iml, jml, lml, tml, 1, 1, totsed) |
---|
1018 | !WRITE(numout,*) "Open erosion effgrid_perc",iml, jml, lml, tml,is_root_prc,fid |
---|
1019 | IF (is_root_prc) CALL flinget(fid, 'effgrid_perc', iml, jml, lml, tml, 1, 1, effgrid) |
---|
1020 | !WRITE(numout,*) "Open erosion fid",iml, jml, lml, tml,is_root_prc,fid |
---|
1021 | IF (is_root_prc) CALL flinclo(fid) |
---|
1022 | ! |
---|
1023 | !WRITE(numout,*) 'lat-forcingdata', lalo(1:4,1) |
---|
1024 | !WRITE(numout,*) 'lon-forcingdata', lalo(1:4,2) |
---|
1025 | WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel),lonrel(1:2,1) |
---|
1026 | WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel),latrel(1,1:2) |
---|
1027 | ! WRITE(numout,*) 'erosion_rate : ', totsed(:,:) |
---|
1028 | ! WRITE(numout,*) 'effgrid : ', effgrid(:,:) |
---|
1029 | ! |
---|
1030 | DO ip=1,iml |
---|
1031 | DO jp=1,jml |
---|
1032 | IF ( totsed(ip,jp).GT. min_sechiba ) THEN |
---|
1033 | mask(ip,jp) = 1 |
---|
1034 | ENDIF |
---|
1035 | ! Resolution in longitude |
---|
1036 | ! |
---|
1037 | coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos ) |
---|
1038 | IF ( ip .EQ. 1 ) THEN |
---|
1039 | resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat |
---|
1040 | ELSEIF ( ip .EQ. iml ) THEN |
---|
1041 | resol_lu(ip,jp,1) = ABS( lonrel(ip,jp) - lonrel(ip-1,jp) ) * pi/180. * R_Earth * coslat |
---|
1042 | ELSE |
---|
1043 | resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip-1,jp) )/2. * pi/180. * R_Earth * coslat |
---|
1044 | ENDIF |
---|
1045 | ! |
---|
1046 | ! Resolution in latitude |
---|
1047 | ! |
---|
1048 | IF ( jp .EQ. 1 ) THEN |
---|
1049 | resol_lu(ip,jp,2) = ABS( latrel(ip,jp) - latrel(ip,jp+1) ) * pi/180. * R_Earth |
---|
1050 | ELSEIF ( jp .EQ. jml ) THEN |
---|
1051 | resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp) ) * pi/180. * R_Earth |
---|
1052 | ELSE |
---|
1053 | resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp+1) )/2. * pi/180. * R_Earth |
---|
1054 | ENDIF |
---|
1055 | ! |
---|
1056 | ENDDO |
---|
1057 | ENDDO |
---|
1058 | !WRITE(numout,*) "lon_resolu=",MAXVAL(resol_lu(:,:,1)), "lat_resolu=",MAXVAL(resol_lu(:,:,2)), iml, jml |
---|
1059 | ! |
---|
1060 | res_lon=lonrel(2,1)-lonrel(1,1) |
---|
1061 | res_lat=latrel(1,1)-latrel(1,2) |
---|
1062 | lon0=MINVAL(lonrel)-res_lon/2.0 |
---|
1063 | lat0=MAXVAL(latrel)+res_lat/2.0 |
---|
1064 | !WRITE(numout,*) 'resolution', res_lon, res_lat |
---|
1065 | !WRITE(numout,*) 'Coordinate of WEST-NORTH cornor (row1, col1)', lon0, lat0 |
---|
1066 | gridarea(:)=zero |
---|
1067 | std_totsed(:)=zero |
---|
1068 | effgrid_perc(:)=zero |
---|
1069 | DO ig =1,nbpt |
---|
1070 | ip=floor((lalo(ig,2)-lon0)/res_lon)+1 |
---|
1071 | jp=floor((lat0-lalo(ig,1))/res_lat)+1 |
---|
1072 | IF ((ip.GT.zero).AND.(ip.LE.iml).AND.(jp.GT.zero).AND.(jp.LE.jml)) THEN |
---|
1073 | gridarea(ig)=resol_lu(ip,jp,1)*resol_lu(ip,jp,2) |
---|
1074 | std_totsed(ig)=MAX(totsed(ip,jp),0.0) |
---|
1075 | effgrid_perc(ig)=MAX(effgrid(ip,jp),0.0) |
---|
1076 | ENDIF !((ip.GT.zero).AND.(ip.LE.iml).AND.(jp.GT.zero).AND.(jp.LE.jml)) THEN |
---|
1077 | !WRITE(numout,*) 'Grid_id:', ig,'longitude',lalo(ig,2),'Latitude',lalo(ig,1) |
---|
1078 | !WRITE(numout,*) 'Erosion map reading', ip,jp,gridarea(ig),std_totsed(ig),effgrid_perc(ig) |
---|
1079 | ENDDO !DO ig =1,nbpt |
---|
1080 | is_root_prc=is_root_prc_old |
---|
1081 | ! |
---|
1082 | ! WRITE(numout,*) 'Effective area', gridarea(:) |
---|
1083 | ! WRITE(numout,*) 'Effective grid percent', effgrid_perc(:) |
---|
1084 | ! WRITE(numout,*) 'standrad sediment loss rate(ton)', std_totsed(:) |
---|
1085 | ! |
---|
1086 | ! Variable that will be used in erosion_main(**) should not be deallocated zhc |
---|
1087 | IF (ALLOCATED(latrel)) DEALLOCATE (latrel) |
---|
1088 | IF (ALLOCATED(lonrel)) DEALLOCATE (lonrel) |
---|
1089 | IF (ALLOCATED(mask)) DEALLOCATE (mask) |
---|
1090 | IF (ALLOCATED(resol_lu)) DEALLOCATE (resol_lu) |
---|
1091 | IF (ALLOCATED(totsed)) DEALLOCATE (totsed) |
---|
1092 | IF (ALLOCATED(effgrid)) DEALLOCATE (effgrid) |
---|
1093 | ! |
---|
1094 | END SUBROUTINE erosionmap |
---|
1095 | ! |
---|
1096 | ! |
---|
1097 | ! |
---|
1098 | END MODULE erosion |
---|
1099 | |
---|