1 | PROGRAM forcesoil |
---|
2 | !--------------------------------------------------------------------- |
---|
3 | ! This program reads a forcing file for STOMATE's soil carbon routine |
---|
4 | ! which was created with STOMATE. |
---|
5 | ! It then integrates STOMATE's soil carbon parameterizations |
---|
6 | ! for a given number of years using this forcing. |
---|
7 | !--------------------------------------------------------------------- |
---|
8 | !- IPSL (2006) |
---|
9 | !- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
10 | USE netcdf |
---|
11 | !- |
---|
12 | USE defprec |
---|
13 | USE constantes |
---|
14 | USE constantes_veg |
---|
15 | USE constantes_co2 |
---|
16 | USE stomate_constants |
---|
17 | USE ioipsl |
---|
18 | USE stomate_soilcarbon |
---|
19 | USE parallel |
---|
20 | !- |
---|
21 | IMPLICIT NONE |
---|
22 | !- |
---|
23 | CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out,var_name |
---|
24 | INTEGER(i_std) :: iim,jjm |
---|
25 | INTEGER(i_std),PARAMETER :: llm = 1 |
---|
26 | INTEGER(i_std) :: kjpindex |
---|
27 | INTEGER(i_std) :: itau_dep,itau_len |
---|
28 | CHARACTER(LEN=30) :: time_str |
---|
29 | INTEGER(i_std) :: ier,iret |
---|
30 | REAL(r_std) :: dt_files |
---|
31 | REAL(r_std) :: date0 |
---|
32 | INTEGER(i_std) :: rest_id_sto |
---|
33 | INTEGER(i_std) :: ncfid |
---|
34 | REAL(r_std) :: dt_force,dt_forcesoil |
---|
35 | INTEGER :: nparan |
---|
36 | INTEGER,PARAMETER :: nparanmax=36 |
---|
37 | REAL(r_std) :: xbid1,xbid2 |
---|
38 | INTEGER(i_std) :: ibid |
---|
39 | INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices |
---|
40 | REAL(r_std),DIMENSION(llm) :: lev |
---|
41 | REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input |
---|
42 | REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: & |
---|
43 | & carbon,control_moist,control_temp |
---|
44 | REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: & |
---|
45 | & lon,lat,resp_hetero_soil,var_3d |
---|
46 | REAL(r_std),DIMENSION(:),ALLOCATABLE :: & |
---|
47 | & x_indices |
---|
48 | REAL(r_std) :: time |
---|
49 | INTEGER :: i,j,m,iatt,iv |
---|
50 | CHARACTER(LEN=400) :: taboo_vars |
---|
51 | REAL(r_std),DIMENSION(1) :: xtmp |
---|
52 | INTEGER(i_std),PARAMETER :: nbvarmax=300 |
---|
53 | INTEGER(i_std) :: nbvar |
---|
54 | CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames |
---|
55 | INTEGER(i_std) :: varnbdim |
---|
56 | INTEGER(i_std),PARAMETER :: varnbdim_max=20 |
---|
57 | INTEGER,DIMENSION(varnbdim_max) :: vardims |
---|
58 | LOGICAL :: l1d |
---|
59 | REAL(r_std) :: x_tmp |
---|
60 | ! clay fraction |
---|
61 | REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay |
---|
62 | !- |
---|
63 | ! string suffix indicating an index |
---|
64 | CHARACTER(LEN=10) :: part_str |
---|
65 | ! |
---|
66 | CHARACTER(LEN=100) :: Cforcing_name |
---|
67 | INTEGER :: Cforcing_id |
---|
68 | INTEGER :: v_id |
---|
69 | |
---|
70 | REAL(r_std),ALLOCATABLE :: clay_loc(:) |
---|
71 | REAL(r_std),ALLOCATABLE :: soilcarbon_input_loc(:,:,:,:) |
---|
72 | REAL(r_std),ALLOCATABLE :: control_temp_loc(:,:,:) |
---|
73 | REAL(r_std),ALLOCATABLE :: control_moist_loc(:,:,:) |
---|
74 | REAL(r_std),ALLOCATABLE :: carbon_loc(:,:,:) |
---|
75 | INTEGER :: ierr |
---|
76 | |
---|
77 | CALL Init_para(.FALSE.) |
---|
78 | |
---|
79 | !- |
---|
80 | ! Stomate's restart files |
---|
81 | !- |
---|
82 | IF (is_root_prc) THEN |
---|
83 | sto_restname_in = 'stomate_start.nc' |
---|
84 | CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) |
---|
85 | WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) |
---|
86 | sto_restname_out = 'stomate_restart.nc' |
---|
87 | CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) |
---|
88 | WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) |
---|
89 | !- |
---|
90 | ! We need to know iim, jjm. |
---|
91 | ! Get them from the restart files themselves. |
---|
92 | !- |
---|
93 | iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, ncfid) |
---|
94 | iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim) |
---|
95 | iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm) |
---|
96 | iret = NF90_CLOSE (ncfid) |
---|
97 | !- |
---|
98 | ! Allocate longitudes and latitudes |
---|
99 | !- |
---|
100 | ALLOCATE (lon(iim,jjm)) |
---|
101 | ALLOCATE (lat(iim,jjm)) |
---|
102 | lon(:,:) = 0.0 |
---|
103 | lat(:,:) = 0.0 |
---|
104 | lev(1) = 0.0 |
---|
105 | !- |
---|
106 | CALL restini & |
---|
107 | & (sto_restname_in, iim, jjm, lon, lat, llm, lev, & |
---|
108 | & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) |
---|
109 | ENDIF |
---|
110 | !- |
---|
111 | ! calendar |
---|
112 | !- |
---|
113 | CALL bcast(date0) |
---|
114 | !!! MM : Ã revoir : choix du calendrier dans forcesoil ?? Il est dans le restart de stomate ! |
---|
115 | ! CALL ioconf_calendar ('noleap') |
---|
116 | CALL ioget_calendar (one_year,one_day) |
---|
117 | |
---|
118 | CALL ioconf_startdate(date0) |
---|
119 | |
---|
120 | IF (is_root_prc) THEN |
---|
121 | !- |
---|
122 | ! open FORCESOIL's forcing file to read some basic info |
---|
123 | !- |
---|
124 | Cforcing_name = 'stomate_Cforcing.nc' |
---|
125 | CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name) |
---|
126 | !- |
---|
127 | ier = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) |
---|
128 | !- |
---|
129 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp) |
---|
130 | kjpindex = NINT(x_tmp) |
---|
131 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp) |
---|
132 | nparan = NINT(x_tmp) |
---|
133 | !- |
---|
134 | ALLOCATE (indices(kjpindex)) |
---|
135 | ALLOCATE (clay(kjpindex)) |
---|
136 | !- |
---|
137 | ALLOCATE (x_indices(kjpindex),stat=ier) |
---|
138 | ier = NF90_INQ_VARID (Cforcing_id,'index',v_id) |
---|
139 | ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices) |
---|
140 | indices(:) = NINT(x_indices(:)) |
---|
141 | DEALLOCATE (x_indices) |
---|
142 | !- |
---|
143 | ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id) |
---|
144 | ier = NF90_GET_VAR (Cforcing_id,v_id,clay) |
---|
145 | !- |
---|
146 | ! time step of forcesoil |
---|
147 | !- |
---|
148 | dt_forcesoil = one_year / FLOAT(nparan) |
---|
149 | WRITE(*,*) 'time step (d): ',dt_forcesoil |
---|
150 | !- |
---|
151 | ! read (and partially write) the restart file |
---|
152 | !- |
---|
153 | ! read and write the variables we do not need |
---|
154 | !- |
---|
155 | taboo_vars = '$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// & |
---|
156 | & '$day_counter$ $dt_days$ $date$ '// & |
---|
157 | & '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// & |
---|
158 | & '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// & |
---|
159 | & '$carbon_11$ $carbon_12$ $carbon_13$' |
---|
160 | !- |
---|
161 | CALL ioget_vname(rest_id_sto, nbvar, varnames) |
---|
162 | !- |
---|
163 | ! read and write some special variables (1D or variables that we need) |
---|
164 | !- |
---|
165 | var_name = 'day_counter' |
---|
166 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
167 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
168 | !- |
---|
169 | var_name = 'dt_days' |
---|
170 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
171 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
172 | !- |
---|
173 | var_name = 'date' |
---|
174 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
175 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
176 | !- |
---|
177 | DO iv=1,nbvar |
---|
178 | !-- check if the variable is to be written here |
---|
179 | IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN |
---|
180 | !---- get variable dimensions, especially 3rd dimension |
---|
181 | CALL ioget_vdim & |
---|
182 | & (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims) |
---|
183 | l1d = ALL(vardims(1:varnbdim) == 1) |
---|
184 | !---- |
---|
185 | ALLOCATE( var_3d(kjpindex,vardims(3)), stat=ier) |
---|
186 | IF (ier /= 0) STOP 'ALLOCATION PROBLEM' |
---|
187 | !---- read it |
---|
188 | IF (l1d) THEN |
---|
189 | CALL restget & |
---|
190 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
191 | & 1, itau_dep, .TRUE., xtmp) |
---|
192 | ELSE |
---|
193 | CALL restget & |
---|
194 | & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), & |
---|
195 | & 1, itau_dep, .TRUE., var_3d, "gather", kjpindex, indices) |
---|
196 | ENDIF |
---|
197 | !---- write it |
---|
198 | IF (l1d) THEN |
---|
199 | CALL restput & |
---|
200 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
201 | & 1, itau_dep, xtmp) |
---|
202 | ELSE |
---|
203 | CALL restput & |
---|
204 | & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), & |
---|
205 | & 1, itau_dep, var_3d, 'scatter', kjpindex, indices) |
---|
206 | ENDIF |
---|
207 | !---- |
---|
208 | DEALLOCATE(var_3d) |
---|
209 | ENDIF |
---|
210 | ENDDO |
---|
211 | !- |
---|
212 | ! read soil carbon |
---|
213 | !- |
---|
214 | ALLOCATE(carbon(kjpindex,ncarb,nvm)) |
---|
215 | carbon(:,:,:) = val_exp |
---|
216 | DO m = 1, nvm |
---|
217 | WRITE (part_str, '(I2)') m |
---|
218 | IF (m<10) part_str(1:1)='0' |
---|
219 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
220 | CALL restget & |
---|
221 | & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, & |
---|
222 | & .TRUE., carbon(:,:,m), 'gather', kjpindex, indices) |
---|
223 | IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = zero |
---|
224 | !-- do not write this variable: it will be modified. |
---|
225 | ENDDO |
---|
226 | !- |
---|
227 | ! Length of run |
---|
228 | !- |
---|
229 | WRITE(time_str,'(a)') '10000Y' |
---|
230 | CALL getin('TIME_LENGTH', time_str) |
---|
231 | ! transform into itau |
---|
232 | CALL tlen2itau(time_str, dt_forcesoil*one_year, date0, itau_len) |
---|
233 | write(*,*) 'Number of time steps to do: ',itau_len |
---|
234 | !- |
---|
235 | ! read the rest of the forcing file and store forcing in an array. |
---|
236 | ! We read an average year. |
---|
237 | !- |
---|
238 | ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan)) |
---|
239 | ALLOCATE(control_temp(kjpindex,nlevs,nparan)) |
---|
240 | ALLOCATE(control_moist(kjpindex,nlevs,nparan)) |
---|
241 | !- |
---|
242 | ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id) |
---|
243 | ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input) |
---|
244 | ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id) |
---|
245 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist) |
---|
246 | ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id) |
---|
247 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp) |
---|
248 | !- |
---|
249 | ier = NF90_CLOSE (Cforcing_id) |
---|
250 | !- |
---|
251 | !MM Problem here with dpu which depends on soil type |
---|
252 | DO iv = 1, nbdl-1 |
---|
253 | ! first 2.0 is dpu |
---|
254 | ! second 2.0 is average |
---|
255 | diaglev(iv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(iv-1) -1) + ( 2**(iv) -1) ) / 2.0 |
---|
256 | ENDDO |
---|
257 | diaglev(nbdl) = 2.0 |
---|
258 | !- |
---|
259 | ! For sequential use only, we must initialize data_para : |
---|
260 | ! |
---|
261 | ! |
---|
262 | ENDIF |
---|
263 | |
---|
264 | CALL bcast(iim) |
---|
265 | CALL bcast(jjm) |
---|
266 | call bcast(kjpindex) |
---|
267 | CALL init_data_para(iim,jjm,kjpindex,indices) |
---|
268 | |
---|
269 | !- |
---|
270 | !- |
---|
271 | ! there we go: time loop |
---|
272 | !- |
---|
273 | CALL bcast(nparan) |
---|
274 | ALLOCATE(clay_loc(nbp_loc)) |
---|
275 | ALLOCATE(soilcarbon_input_loc(nbp_loc,ncarb,nvm,nparan)) |
---|
276 | ALLOCATE(control_temp_loc(nbp_loc,nlevs,nparan)) |
---|
277 | ALLOCATE(control_moist_loc(nbp_loc,nlevs,nparan)) |
---|
278 | ALLOCATE(carbon_loc(nbp_loc,ncarb,nvm)) |
---|
279 | ALLOCATE(resp_hetero_soil(nbp_loc,nvm)) |
---|
280 | iatt = 0 |
---|
281 | |
---|
282 | CALL bcast(itau_len) |
---|
283 | CALL bcast(nparan) |
---|
284 | CALL bcast(dt_forcesoil) |
---|
285 | CALL Scatter(clay,clay_loc) |
---|
286 | CALL Scatter(soilcarbon_input,soilcarbon_input_loc) |
---|
287 | CALL Scatter(control_temp,control_temp_loc) |
---|
288 | CALL Scatter(control_moist,control_moist_loc) |
---|
289 | CALL Scatter(carbon,carbon_loc) |
---|
290 | |
---|
291 | DO i=1,itau_len |
---|
292 | iatt = iatt+1 |
---|
293 | IF (iatt > nparan) iatt = 1 |
---|
294 | CALL soilcarbon & |
---|
295 | & (nbp_loc, dt_forcesoil, clay_loc, & |
---|
296 | & soilcarbon_input_loc(:,:,:,iatt), & |
---|
297 | & control_temp_loc(:,:,iatt), control_moist_loc(:,:,iatt), & |
---|
298 | & carbon_loc, resp_hetero_soil) |
---|
299 | ENDDO |
---|
300 | |
---|
301 | CALL Gather(carbon_loc,carbon) |
---|
302 | !- |
---|
303 | ! write new carbon into restart file |
---|
304 | !- |
---|
305 | IF (is_root_prc) THEN |
---|
306 | DO m=1,nvm |
---|
307 | WRITE (part_str, '(I2)') m |
---|
308 | IF (m<10) part_str(1:1)='0' |
---|
309 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
310 | CALL restput & |
---|
311 | & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, & |
---|
312 | & carbon(:,:,m), 'scatter', kjpindex, indices) |
---|
313 | ENDDO |
---|
314 | !- |
---|
315 | CALL getin_dump |
---|
316 | CALL restclo |
---|
317 | ENDIF |
---|
318 | #ifdef CPP_PARA |
---|
319 | CALL MPI_FINALIZE(ierr) |
---|
320 | #endif |
---|
321 | !-------------------- |
---|
322 | END PROGRAM forcesoil |
---|