| 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 | |
|---|
