1 | ! |
---|
2 | ! Daily update of leaf area index etc. |
---|
3 | ! |
---|
4 | !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/slowproc.f90,v 1.48 2010/04/20 14:12:04 ssipsl Exp $ |
---|
5 | !! IPSL (2006) |
---|
6 | !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
7 | ! |
---|
8 | ! |
---|
9 | MODULE slowproc |
---|
10 | |
---|
11 | ! modules used: |
---|
12 | |
---|
13 | USE defprec |
---|
14 | USE constantes |
---|
15 | USE constantes_veg |
---|
16 | USE constantes_co2 |
---|
17 | USE ioipsl |
---|
18 | USE sechiba_io |
---|
19 | USE interpol_help |
---|
20 | USE stomate |
---|
21 | USE stomate_constants |
---|
22 | USE grid |
---|
23 | USE parallel |
---|
24 | ! USE Write_Field_p |
---|
25 | |
---|
26 | IMPLICIT NONE |
---|
27 | |
---|
28 | PRIVATE |
---|
29 | PUBLIC slowproc_main,slowproc_clear |
---|
30 | |
---|
31 | ! To use OLD or NEW iterpollation schemes : |
---|
32 | INTERFACE slowproc_interlai |
---|
33 | MODULE PROCEDURE slowproc_interlai_OLD, slowproc_interlai_NEW |
---|
34 | END INTERFACE |
---|
35 | INTERFACE slowproc_interpol |
---|
36 | MODULE PROCEDURE slowproc_interpol_OLD, slowproc_interpol_NEW |
---|
37 | END INTERFACE |
---|
38 | INTERFACE slowproc_interpol_g |
---|
39 | MODULE PROCEDURE slowproc_interpol_OLD_g, slowproc_interpol_NEW_g |
---|
40 | END INTERFACE |
---|
41 | |
---|
42 | ! |
---|
43 | ! variables used inside slowproc module : declaration and initialisation |
---|
44 | ! |
---|
45 | |
---|
46 | LOGICAL, SAVE :: l_first_slowproc = .TRUE.!! Initialisation has to be done one time |
---|
47 | REAL(r_std), SAVE :: dt_slow !! time step of slow processes and STOMATE |
---|
48 | ! |
---|
49 | REAL(r_std), SAVE :: clayfraction_default = 0.2 |
---|
50 | |
---|
51 | INTEGER(i_std) , SAVE :: veget_update=0 !! update frequency in years for landuse |
---|
52 | INTEGER(i_std) , SAVE :: veget_year_orig=0 !! first year for landuse |
---|
53 | LOGICAL, SAVE :: land_use = .FALSE. ! Land Use |
---|
54 | LOGICAL, SAVE :: veget_reinit=.FALSE. !! To change LAND USE file in a run. |
---|
55 | ! |
---|
56 | LOGICAL, SAVE :: read_lai = .FALSE. ! Lai Map |
---|
57 | LOGICAL, SAVE :: old_lai = .FALSE. ! Old Lai Map interpolation |
---|
58 | LOGICAL, SAVE :: impveg = .FALSE. |
---|
59 | LOGICAL, SAVE :: old_veget = .FALSE. ! Old veget Map interpolation |
---|
60 | ! |
---|
61 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: clayfraction |
---|
62 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laimap |
---|
63 | ! |
---|
64 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_nextyear !! next year fraction of vegetation type |
---|
65 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_nextyear !! next year fraction of ice+lakes+cities+... |
---|
66 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio_nextyear !! next year total fraction of ice+lakes+cities+... |
---|
67 | |
---|
68 | CONTAINS |
---|
69 | |
---|
70 | SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, & |
---|
71 | ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & |
---|
72 | IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltype, & |
---|
73 | t2m, t2m_min, temp_sol, stempdiag, & |
---|
74 | humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & |
---|
75 | deadleaf_cover, & |
---|
76 | assim_param, & |
---|
77 | lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, & |
---|
78 | rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & |
---|
79 | co2_flux) |
---|
80 | |
---|
81 | |
---|
82 | ! interface description |
---|
83 | ! input scalar |
---|
84 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number |
---|
85 | INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid |
---|
86 | INTEGER(i_std),INTENT (in) :: kjpindex !! Domain size |
---|
87 | REAL(r_std),INTENT (in) :: dtradia !! Time step |
---|
88 | REAL(r_std),INTENT (in) :: date0 !! Initial date |
---|
89 | LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read |
---|
90 | LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write |
---|
91 | LOGICAL, INTENT(in) :: ldforcing_write !! Logical for _forcing_ file to write |
---|
92 | LOGICAL, INTENT(in) :: ldcarbon_write !! Logical for _carbon_forcing_ file to write |
---|
93 | INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier |
---|
94 | INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier |
---|
95 | INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file file identifier |
---|
96 | INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file file identifier |
---|
97 | INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier |
---|
98 | ! input fields |
---|
99 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map |
---|
100 | INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map |
---|
101 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) |
---|
102 | INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land |
---|
103 | REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) |
---|
104 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid |
---|
105 | REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity ("moisture stress") |
---|
106 | REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m !! 2 m air temperature (K) |
---|
107 | REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m_min !! min. 2 m air temp. during forcing time step (K) |
---|
108 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature |
---|
109 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature |
---|
110 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil moisture |
---|
111 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: litterhumdiag !! Litter humidity |
---|
112 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation |
---|
113 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation |
---|
114 | REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: gpp !! GPP (gC/(m**2 of total ground)/time step) |
---|
115 | ! output scalar |
---|
116 | ! output fields |
---|
117 | REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux in gC/m**2 of average ground/second |
---|
118 | ! modified scalar |
---|
119 | ! modified fields |
---|
120 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: lai !! Surface foliaire |
---|
121 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: height !! height (m) |
---|
122 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget !! Fraction of vegetation type |
---|
123 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout):: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh |
---|
124 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget_max !! Max fraction of vegetation type |
---|
125 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh |
---|
126 | REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout):: soiltype !! fraction of soil type |
---|
127 | REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (inout):: assim_param!! min+max+opt temps, vcmax, vjmax for photosynthesis |
---|
128 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: deadleaf_cover !! fraction of soil covered by dead leaves |
---|
129 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: qsintmax !! Maximum water on vegetation for interception |
---|
130 | |
---|
131 | ! local declaration |
---|
132 | INTEGER(i_std) :: j, jv |
---|
133 | INTEGER(i_std), SAVE :: lcanop !! soil level used for LAI |
---|
134 | REAL(r_std) :: tmp_day(1) |
---|
135 | CHARACTER(LEN=80) :: var_name !! To store variables names for I/O |
---|
136 | |
---|
137 | ! |
---|
138 | REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_maint !! autotrophic resp. (gC/(m**2 of total ground)/time step) |
---|
139 | REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_hetero !! heterotrophic resp. (gC/(m**2 of total ground)/time step) |
---|
140 | REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_growth !! growth resp. (gC/(m**2 of total ground)/time step) |
---|
141 | REAL(r_std), DIMENSION(kjpindex,nvm) :: npp !! Net Ecosystem Exchange (gC/(m**2 of total ground)/time step) |
---|
142 | ! |
---|
143 | INTEGER(i_std) , SAVE :: veget_year !! year for landuse |
---|
144 | LOGICAL :: land_use_updated |
---|
145 | ! |
---|
146 | LOGICAL, PARAMETER :: check = .FALSE. |
---|
147 | |
---|
148 | REAL(r_std), SAVE :: sec_old = 0. |
---|
149 | ! |
---|
150 | ! do initialisation |
---|
151 | ! |
---|
152 | IF (l_first_slowproc) THEN |
---|
153 | |
---|
154 | ! |
---|
155 | ! 1.1 allocation, file definitions. Restart file read for Sechiba. Set flags. |
---|
156 | ! |
---|
157 | |
---|
158 | IF (long_print) WRITE (numout,*) ' l_first_slowproc : call slowproc_init ' |
---|
159 | CALL slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & |
---|
160 | & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,& |
---|
161 | & veget_update, veget_year) |
---|
162 | ! |
---|
163 | ! Time step in days for stomate |
---|
164 | dt_days = dt_slow / one_day |
---|
165 | ! |
---|
166 | resp_maint(:,:) = zero |
---|
167 | resp_hetero(:,:) = zero |
---|
168 | resp_growth(:,:) = zero |
---|
169 | ! |
---|
170 | ! 1.2 check time step |
---|
171 | ! |
---|
172 | IF ( dt_slow .LT. dtradia ) THEN |
---|
173 | WRITE(numout,*) 'slow_processes: time step smaller than forcing time step.' |
---|
174 | STOP 'slowproc_main' |
---|
175 | ENDIF |
---|
176 | ! |
---|
177 | IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN |
---|
178 | ! |
---|
179 | ! 1.3 call STOMATE for initialization |
---|
180 | ! |
---|
181 | CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & |
---|
182 | ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & |
---|
183 | IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & |
---|
184 | t2m, t2m_min, temp_sol, stempdiag, & |
---|
185 | humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & |
---|
186 | deadleaf_cover, & |
---|
187 | assim_param, & |
---|
188 | lai, height, veget, veget_max, qsintmax, & |
---|
189 | veget_nextyear, totfrac_nobio_nextyear, & |
---|
190 | hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & |
---|
191 | co2_flux,resp_maint,resp_hetero,resp_growth) |
---|
192 | ! |
---|
193 | ENDIF |
---|
194 | ! |
---|
195 | IF ( .NOT. control%ok_stomate ) THEN |
---|
196 | ! |
---|
197 | ! 1.4 initialize some variables |
---|
198 | ! STOMATE diagnoses some variables for SECHIBA: height, deadleaf_cover, etc. |
---|
199 | ! IF SECHIBA is not coupled to STOMATE, then we must set these values otherwise. |
---|
200 | ! |
---|
201 | CALL slowproc_derivvar (kjpindex, veget, lai, & |
---|
202 | qsintmax, deadleaf_cover, assim_param, height) |
---|
203 | ENDIF |
---|
204 | |
---|
205 | RETURN |
---|
206 | |
---|
207 | ENDIF |
---|
208 | ! |
---|
209 | !!$ ! Land USE for next year |
---|
210 | !!$ land_use_updated=.FALSE. |
---|
211 | !!$ IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN |
---|
212 | !!$ ! if next iteration is divisibled by veget_update |
---|
213 | !!$ IF ( mod(kjit+1, veget_update) .le. min_sechiba) THEN |
---|
214 | !!$ ! |
---|
215 | !!$ veget_year=veget_year+veget_year_add |
---|
216 | !!$ WRITE(numout,*) 'We are updating veget for year =' , veget_year |
---|
217 | !!$ ! |
---|
218 | !!$ ! Save veget |
---|
219 | !!$ ! |
---|
220 | !!$ veget_lastyear(:,:) = veget_max(:,:) |
---|
221 | !!$ ! |
---|
222 | !!$ CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & |
---|
223 | !!$ & veget_max, frac_nobio, veget_year) |
---|
224 | !!$ |
---|
225 | !!$ land_use_updated=.TRUE. |
---|
226 | !!$ ENDIF |
---|
227 | !!$ ENDIF |
---|
228 | ! |
---|
229 | ! prepares restart file for the next simulation |
---|
230 | ! |
---|
231 | IF (ldrestart_write) THEN |
---|
232 | |
---|
233 | IF (long_print) WRITE (numout,*) ' we have to complete restart file with SLOWPROC variables ' |
---|
234 | |
---|
235 | tmp_day(1) = day_counter |
---|
236 | var_name= 'day_counter' |
---|
237 | IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day) |
---|
238 | |
---|
239 | var_name= 'veget' |
---|
240 | CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget, 'scatter', nbp_glo, index_g) |
---|
241 | ! |
---|
242 | var_name= 'veget_max' |
---|
243 | CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget_max, 'scatter', nbp_glo, index_g) |
---|
244 | ! |
---|
245 | var_name= 'lai' |
---|
246 | CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, lai, 'scatter', nbp_glo, index_g) |
---|
247 | ! |
---|
248 | var_name= 'frac_nobio' |
---|
249 | CALL restput_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter', nbp_glo, index_g) |
---|
250 | ! |
---|
251 | var_name= 'soiltype_frac' |
---|
252 | CALL restput_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, soiltype, 'scatter', nbp_glo, index_g) |
---|
253 | ! |
---|
254 | var_name= 'clay_frac' |
---|
255 | CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, clayfraction, 'scatter', nbp_glo, index_g) |
---|
256 | ! |
---|
257 | ! The height of the vegetation could in principle be recalculated at the beginning of the run. |
---|
258 | ! However, this is very tedious, as many special cases have to be taken into account. This variable |
---|
259 | ! is therefore saved in the restart file. |
---|
260 | var_name= 'height' |
---|
261 | CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, height, 'scatter', nbp_glo, index_g) |
---|
262 | ! |
---|
263 | IF (read_lai) THEN |
---|
264 | var_name= 'laimap' |
---|
265 | CALL restput_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, laimap) |
---|
266 | ENDIF |
---|
267 | ! |
---|
268 | IF (land_use) THEN |
---|
269 | tmp_day(1) = REAL(veget_year,r_std) |
---|
270 | var_name='veget_year' |
---|
271 | IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day) |
---|
272 | ENDIF |
---|
273 | ! |
---|
274 | ! call STOMATE to write RESTART files |
---|
275 | ! |
---|
276 | IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN |
---|
277 | CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & |
---|
278 | ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & |
---|
279 | IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & |
---|
280 | t2m, t2m_min, temp_sol, stempdiag, & |
---|
281 | humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & |
---|
282 | deadleaf_cover, & |
---|
283 | assim_param, & |
---|
284 | lai, height, veget, veget_max, qsintmax, & |
---|
285 | veget_nextyear, totfrac_nobio_nextyear, & |
---|
286 | hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & |
---|
287 | co2_flux,resp_maint,resp_hetero,resp_growth) |
---|
288 | ENDIF |
---|
289 | |
---|
290 | RETURN |
---|
291 | |
---|
292 | END IF |
---|
293 | ! |
---|
294 | ! update day counter |
---|
295 | day_counter = day_counter + dtradia |
---|
296 | |
---|
297 | IF (check) WRITE(*,*) "slowproc: day_counter 3",day_counter |
---|
298 | ! |
---|
299 | ! |
---|
300 | ! 1. Compute date |
---|
301 | ! |
---|
302 | ! Test each day and assert all slow processes (days and years) |
---|
303 | ! |
---|
304 | IF ( sec_old >= one_day - dtradia .AND. sec >= 0. ) THEN |
---|
305 | ! |
---|
306 | ! reset counter |
---|
307 | day_counter = zero |
---|
308 | IF (check) WRITE(*,*) "slowproc: day_counter 2",day_counter |
---|
309 | ! |
---|
310 | ! Active slow processes |
---|
311 | do_slow = .TRUE. |
---|
312 | ! |
---|
313 | ! count days |
---|
314 | date = date + nint(dt_days) |
---|
315 | IF (check) WRITE(numout,*) "New date : ",date, 'year_length ',year_length,kjit |
---|
316 | ! |
---|
317 | ! is one year over? |
---|
318 | ! EndOfYear must be true once per year |
---|
319 | ! during a call of stomate_season. |
---|
320 | ! |
---|
321 | IF ( month == 1 .AND. day == 1 .AND. sec .LT. dtradia ) THEN |
---|
322 | EndOfYear = .TRUE. |
---|
323 | ELSE |
---|
324 | EndOfYear = .FALSE. |
---|
325 | ENDIF |
---|
326 | |
---|
327 | ELSE |
---|
328 | do_slow = .FALSE. |
---|
329 | EndOfYear = .FALSE. |
---|
330 | ENDIF |
---|
331 | sec_old = sec |
---|
332 | |
---|
333 | IF ( EndOfYear ) THEN |
---|
334 | WRITE(numout,*) 'slowproc: EndOfYear is activated.' |
---|
335 | ENDIF |
---|
336 | |
---|
337 | ! Land USE for next year |
---|
338 | land_use_updated=.FALSE. |
---|
339 | IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN |
---|
340 | ! if next iteration is divisibled by veget_update |
---|
341 | IF ( EndOfYear ) THEN |
---|
342 | ! |
---|
343 | veget_year = veget_year + 1 |
---|
344 | ! |
---|
345 | IF ( MOD(veget_year - veget_year_orig, veget_update) == 0 ) THEN |
---|
346 | |
---|
347 | WRITE(numout,*) 'We are updating land use veget for year =' , veget_year |
---|
348 | ! |
---|
349 | CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & |
---|
350 | & veget_max, frac_nobio, veget_nextyear, frac_nobio_nextyear, veget_year) |
---|
351 | ! |
---|
352 | ! If some PFT has disapeared in new map |
---|
353 | WHERE(veget_nextyear(:,:).LT.veget(:,:)) |
---|
354 | veget(:,:) = veget_nextyear(:,:) |
---|
355 | ENDWHERE |
---|
356 | ! |
---|
357 | DO j = 1, nnobio |
---|
358 | totfrac_nobio_nextyear(:) = totfrac_nobio_nextyear(:) + frac_nobio_nextyear(:,j) |
---|
359 | ENDDO |
---|
360 | land_use_updated=.TRUE. |
---|
361 | ! |
---|
362 | ENDIF |
---|
363 | ! |
---|
364 | ENDIF |
---|
365 | ENDIF ! Land Use part |
---|
366 | IF (EndOfYear .AND. .NOT. land_use_updated) THEN |
---|
367 | lcchange=.FALSE. |
---|
368 | ENDIF |
---|
369 | |
---|
370 | IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN |
---|
371 | ! |
---|
372 | ! 2. call STOMATE, either because we want to keep track of long-term variables or |
---|
373 | ! because STOMATE is activated |
---|
374 | ! |
---|
375 | CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & |
---|
376 | ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & |
---|
377 | IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & |
---|
378 | t2m, t2m_min, temp_sol, stempdiag, & |
---|
379 | humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & |
---|
380 | deadleaf_cover, & |
---|
381 | assim_param, & |
---|
382 | lai, height, veget, veget_max, qsintmax, & |
---|
383 | veget_nextyear, totfrac_nobio_nextyear, & |
---|
384 | hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & |
---|
385 | co2_flux,resp_maint,resp_hetero,resp_growth) |
---|
386 | IF ( control%ok_stomate .AND. control%ok_sechiba ) THEN |
---|
387 | CALL histwrite(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg) |
---|
388 | CALL histwrite(hist_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg) |
---|
389 | CALL histwrite(hist_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg) |
---|
390 | npp(:,1)=zero |
---|
391 | DO j = 2,nvm |
---|
392 | npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j) |
---|
393 | ENDDO |
---|
394 | CALL histwrite(hist_id, 'npp', kjit, npp, kjpindex*nvm, indexveg) |
---|
395 | ENDIF |
---|
396 | IF ( hist2_id > 0 ) THEN |
---|
397 | IF ( control%ok_stomate ) THEN |
---|
398 | CALL histwrite(hist2_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg) |
---|
399 | CALL histwrite(hist2_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg) |
---|
400 | CALL histwrite(hist2_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg) |
---|
401 | CALL histwrite(hist2_id, 'npp', kjit, npp, kjpindex*nvm, indexveg) |
---|
402 | ENDIF |
---|
403 | ENDIF |
---|
404 | ENDIF |
---|
405 | ! |
---|
406 | IF ( .NOT. control%ok_stomate ) THEN |
---|
407 | |
---|
408 | ! |
---|
409 | ! 2 STOMATE is not activated: we have to guess lai etc. |
---|
410 | ! |
---|
411 | |
---|
412 | ! |
---|
413 | ! 2.2 test if we have work to do |
---|
414 | ! |
---|
415 | |
---|
416 | IF ( do_slow .OR. land_use_updated ) THEN |
---|
417 | ! |
---|
418 | ! 2.2.1 do daily processes if necessary |
---|
419 | ! |
---|
420 | IF (long_print) WRITE (numout,*) 'slowproc_main : We update the daily variables' |
---|
421 | |
---|
422 | ! 2.2.2 updates lai |
---|
423 | CALL slowproc_lai (kjpindex, lcanop,stempdiag, & |
---|
424 | lalo,resolution,lai,month,day,read_lai,laimap) |
---|
425 | ! |
---|
426 | IF (land_use_updated) THEN |
---|
427 | veget_max(:,:)=veget_nextyear(:,:) |
---|
428 | frac_nobio(:,:)=frac_nobio_nextyear(:,:) |
---|
429 | ENDIF |
---|
430 | ! |
---|
431 | ! 2.2.3 updates veget |
---|
432 | CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) |
---|
433 | totfrac_nobio(:) = zero |
---|
434 | DO jv = 1, nnobio |
---|
435 | totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv) |
---|
436 | ENDDO |
---|
437 | |
---|
438 | ! 2.2.4 updates qsintmax and other derived variables |
---|
439 | CALL slowproc_derivvar (kjpindex, veget, lai, & |
---|
440 | qsintmax, deadleaf_cover, assim_param, height) |
---|
441 | ELSE |
---|
442 | ! |
---|
443 | IF (land_use_updated) THEN |
---|
444 | frac_nobio(:,:)=frac_nobio_nextyear(:,:) |
---|
445 | ENDIF |
---|
446 | ! |
---|
447 | END IF |
---|
448 | |
---|
449 | ! |
---|
450 | ! 2.3 some output fields |
---|
451 | ! |
---|
452 | |
---|
453 | co2_flux(:,:) = zero |
---|
454 | |
---|
455 | ENDIF |
---|
456 | |
---|
457 | IF (long_print) WRITE (numout,*) ' slowproc_main done ' |
---|
458 | |
---|
459 | END SUBROUTINE slowproc_main |
---|
460 | !! |
---|
461 | !! |
---|
462 | !! |
---|
463 | SUBROUTINE slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & |
---|
464 | & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,& |
---|
465 | & veget_update, veget_year) |
---|
466 | |
---|
467 | ! interface description |
---|
468 | ! input scalar |
---|
469 | INTEGER(i_std), INTENT (in) :: kjit !! Time step number |
---|
470 | LOGICAL, INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read |
---|
471 | REAL(r_std),INTENT (in) :: dtradia !! Time step |
---|
472 | REAL(r_std), INTENT (in) :: date0 !! intial date |
---|
473 | INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size |
---|
474 | INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier |
---|
475 | ! input fields |
---|
476 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map |
---|
477 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) |
---|
478 | INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land |
---|
479 | REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) |
---|
480 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid |
---|
481 | ! output scalar |
---|
482 | INTEGER(i_std), INTENT(out) :: lcanop !! soil level used for LAI |
---|
483 | INTEGER(i_std), INTENT(out) :: veget_update !! update frequency in timesteps for landuse |
---|
484 | INTEGER(i_std), INTENT(out) :: veget_year !! first year for landuse |
---|
485 | ! output fields |
---|
486 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Surface foliere |
---|
487 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type |
---|
488 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice,lakes,cities, ... |
---|
489 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... |
---|
490 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Max fraction of vegetation type |
---|
491 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! Height of vegetation |
---|
492 | REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltype !! fraction of soil type |
---|
493 | ! local declaration |
---|
494 | REAL(r_std) :: tmp_day(1) |
---|
495 | REAL(r_std) :: tmp_veget_year(1) |
---|
496 | REAL(r_std) :: zcanop !! soil depth taken for canopy |
---|
497 | INTEGER(i_std) :: vtmp(1) |
---|
498 | REAL(r_std), DIMENSION(nbdl) :: zsoil !! soil depths at diagnostic levels |
---|
499 | INTEGER(i_std) :: j,l !! Index |
---|
500 | CHARACTER(LEN=80) :: var_name !! To store variables names for I/O |
---|
501 | INTEGER(i_std) :: ji, jv, ier |
---|
502 | LOGICAL, INTENT(out) :: read_lai |
---|
503 | REAL(r_std) :: frac_nobio1 !! temporary |
---|
504 | REAL(r_std) :: stempdiag_bid !! only needed for an initial LAI |
---|
505 | !! if there is no restart file |
---|
506 | REAL(r_std), DIMENSION(kjpindex,nbdl) :: stempdiag2_bid !! matrix to store stempdiag_bid |
---|
507 | CHARACTER(LEN=4) :: vegsoil_dist !! Flag to choose the soil/vegetation distribution |
---|
508 | ! |
---|
509 | CHARACTER(LEN=30), SAVE :: veget_str !! update frequency for landuse |
---|
510 | REAL(r_std),DIMENSION (nbp_glo,nnobio) :: frac_nobio_g !! Fraction of ice, lakes, cities etc. in the mesh (global) |
---|
511 | REAL(r_std),DIMENSION (nbp_glo,nvm) :: veget_max_g !! Fraction of vegetation type (globa) |
---|
512 | |
---|
513 | LOGICAL, PARAMETER :: check = .FALSE. |
---|
514 | ! |
---|
515 | ! 1 allocation |
---|
516 | ! |
---|
517 | |
---|
518 | ALLOCATE (clayfraction(kjpindex),stat=ier) |
---|
519 | IF (ier.NE.0) THEN |
---|
520 | WRITE (numout,*) ' error in clayfraction allocation. We stop. We need kjpindex words = ',kjpindex |
---|
521 | STOP 'slowproc_init' |
---|
522 | END IF |
---|
523 | clayfraction(:)=undef_sechiba |
---|
524 | ! |
---|
525 | veget_max(:,1) = un |
---|
526 | veget_max(:,2:nvm) = zero |
---|
527 | frac_nobio(:,:) = zero |
---|
528 | totfrac_nobio(:) = zero |
---|
529 | ! |
---|
530 | ier=-1 |
---|
531 | ALLOCATE(veget_nextyear(kjpindex, nvm), STAT=ier) |
---|
532 | IF (ier/=0) THEN |
---|
533 | WRITE(numout,*) "ERROR IN ALLOCATION of veget_nextyear : ",ier |
---|
534 | STOP |
---|
535 | ENDIF |
---|
536 | veget_nextyear(:,1) = un |
---|
537 | veget_nextyear(:,2:nvm) = zero |
---|
538 | ! |
---|
539 | ier=-1 |
---|
540 | ALLOCATE(frac_nobio_nextyear(kjpindex, nnobio), STAT=ier) |
---|
541 | IF (ier/=0) THEN |
---|
542 | PRINT *,"ERROR IN ALLOCATION of frac_nobio_nextyear : ",ier |
---|
543 | STOP |
---|
544 | ENDIF |
---|
545 | frac_nobio_nextyear(:,:) = zero |
---|
546 | ! |
---|
547 | ier=-1 |
---|
548 | ALLOCATE(totfrac_nobio_nextyear(kjpindex), STAT=ier) |
---|
549 | IF (ier/=0) THEN |
---|
550 | PRINT *,"ERROR IN ALLOCATION of totfrac_nobio_nextyear : ",ier |
---|
551 | STOP |
---|
552 | ENDIF |
---|
553 | !MM must be corrected when nnobio > 1 !! |
---|
554 | totfrac_nobio_nextyear(:) = nnobio*un |
---|
555 | ! |
---|
556 | ! 2 read restart file |
---|
557 | ! |
---|
558 | |
---|
559 | var_name= 'day_counter' |
---|
560 | CALL ioconf_setatt('UNITS', 'd') |
---|
561 | CALL ioconf_setatt('LONG_NAME','Fraction of computed day') |
---|
562 | IF (is_root_prc) THEN |
---|
563 | CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_day) |
---|
564 | IF (tmp_day(1) == val_exp) THEN |
---|
565 | day_counter = zero |
---|
566 | ELSE |
---|
567 | day_counter = tmp_day(1) |
---|
568 | ENDIF |
---|
569 | ENDIF |
---|
570 | CALL bcast(day_counter) |
---|
571 | |
---|
572 | ! get restart value if none were found in the restart file |
---|
573 | ! |
---|
574 | !Config Key = SECHIBA_DAY |
---|
575 | !Config Desc = Time within the day simulated |
---|
576 | !Config Def = 0.0 |
---|
577 | !Config Help = This is the time spent simulating the current day. This variable is |
---|
578 | !Config prognostic as it will trigger all the computations which are |
---|
579 | !Config only done once a day. |
---|
580 | ! |
---|
581 | CALL setvar_p (day_counter, val_exp, 'SECHIBA_DAY', 0.0_r_std) |
---|
582 | ! |
---|
583 | !Config Key = LAI_MAP |
---|
584 | !Config Desc = Read the LAI map |
---|
585 | !Config Def = n |
---|
586 | !Config Help = It is possible to read a 12 month LAI map which will |
---|
587 | !Config then be interpolated to daily values as needed. |
---|
588 | ! |
---|
589 | read_lai = .FALSE. |
---|
590 | CALL getin_p('LAI_MAP',read_lai) |
---|
591 | ! |
---|
592 | var_name= 'veget' |
---|
593 | CALL ioconf_setatt('UNITS', '-') |
---|
594 | CALL ioconf_setatt('LONG_NAME','Vegetation fraction') |
---|
595 | CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g) |
---|
596 | ! |
---|
597 | var_name= 'veget_max' |
---|
598 | CALL ioconf_setatt('UNITS', '-') |
---|
599 | CALL ioconf_setatt('LONG_NAME','Maximum vegetation fraction') |
---|
600 | CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g) |
---|
601 | |
---|
602 | ! |
---|
603 | frac_nobio(:,:) = val_exp |
---|
604 | var_name= 'frac_nobio' |
---|
605 | CALL ioconf_setatt('UNITS', '-') |
---|
606 | CALL ioconf_setatt('LONG_NAME','Special soil type fraction') |
---|
607 | CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g) |
---|
608 | ! |
---|
609 | !Config Key = LAND_USE |
---|
610 | !Config Desc = Read a land_use vegetation map |
---|
611 | !Config Def = n |
---|
612 | !Config Help = pft values are needed, max time axis is 293 |
---|
613 | ! |
---|
614 | land_use = .FALSE. |
---|
615 | veget_update=0 |
---|
616 | CALL getin_p('LAND_USE',land_use) |
---|
617 | IF (land_use) THEN |
---|
618 | ! |
---|
619 | !Config Key = VEGET_YEAR |
---|
620 | !Config Desc = Year of the land_use vegetation map to be read (0 == NO TIME AXIS) |
---|
621 | !Config If = LAND_USE |
---|
622 | !Config Def = 282 |
---|
623 | !Config Help = First year for landuse vegetation (2D map by pft). |
---|
624 | !Config Help If VEGET_YEAR == 0, this means there is no time axis. |
---|
625 | ! |
---|
626 | veget_year_orig=282 |
---|
627 | CALL getin_p('VEGET_YEAR', veget_year_orig) |
---|
628 | ! |
---|
629 | !Config Key = VEGET_REINIT |
---|
630 | !Config Desc = booleen to indicate that a new LAND USE file will be used. |
---|
631 | !Config If = LAND_USE |
---|
632 | !Config Def = n |
---|
633 | !Config Help = The parameter is used to bypass veget_year count |
---|
634 | !Config Help and reinitialize it with VEGET_YEAR parameter. |
---|
635 | !Config Help Then it is possible to change LAND USE file. |
---|
636 | ! |
---|
637 | veget_reinit = .FALSE. |
---|
638 | CALL getin_p('VEGET_REINIT', veget_reinit) |
---|
639 | ! |
---|
640 | ! |
---|
641 | var_name= 'veget_year' |
---|
642 | CALL ioconf_setatt('UNITS', '-') |
---|
643 | CALL ioconf_setatt('LONG_NAME','Last year get in Land Use file.') |
---|
644 | IF (is_root_prc) THEN |
---|
645 | CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_veget_year) |
---|
646 | ! |
---|
647 | IF (tmp_veget_year(1) == val_exp) THEN |
---|
648 | veget_year=veget_year_orig |
---|
649 | ELSE |
---|
650 | IF (veget_reinit) THEN |
---|
651 | veget_year=veget_year_orig |
---|
652 | ELSE |
---|
653 | veget_year=INT(tmp_veget_year(1)) |
---|
654 | ENDIF |
---|
655 | ENDIF |
---|
656 | ENDIF |
---|
657 | CALL bcast(veget_year) |
---|
658 | ! |
---|
659 | !Config Key = VEGET_UPDATE |
---|
660 | !Config Desc = Update vegetation frequency |
---|
661 | !Config If = LAND_USE |
---|
662 | !Config Def = 0Y |
---|
663 | !Config Help = The veget datas will be update each this time step. |
---|
664 | ! |
---|
665 | veget_update=0 |
---|
666 | WRITE(veget_str,'(a)') '0Y' |
---|
667 | CALL getin_p('VEGET_UPDATE', veget_str) |
---|
668 | ! |
---|
669 | ! |
---|
670 | l=INDEX(TRIM(veget_str),'Y') |
---|
671 | READ(veget_str(1:(l-1)),"(I2.2)") veget_update |
---|
672 | WRITE(numout,*) "Update frequency for land use in years :",veget_update |
---|
673 | ! |
---|
674 | !Config Key = LAND_COVER_CHANGE |
---|
675 | !Config Desc = treat land use modifications |
---|
676 | !Config If = LAND_USE |
---|
677 | !Config Def = y |
---|
678 | !Config Help = With this variable, you can use a Land Use map |
---|
679 | !Config to simulate anthropic modifications such as |
---|
680 | !Config deforestation. |
---|
681 | ! |
---|
682 | lcchange = .TRUE. |
---|
683 | CALL getin_p('LAND_COVER_CHANGE', lcchange) |
---|
684 | IF ( veget_update == 0 .AND. lcchange ) THEN |
---|
685 | CALL ipslerr (2,'slowproc_init', & |
---|
686 | & 'You have asked for LAND_COVER_CHANGE activated with VEGET_UPDATE = 0Y.',& |
---|
687 | & 'We can''t use this land cover change model if veget is not updated.', & |
---|
688 | & 'We have disabled it.') |
---|
689 | lcchange=.FALSE. |
---|
690 | ENDIF |
---|
691 | |
---|
692 | ENDIF |
---|
693 | ! |
---|
694 | var_name= 'soiltype_frac' |
---|
695 | CALL ioconf_setatt('UNITS', '-') |
---|
696 | CALL ioconf_setatt('LONG_NAME','Fraction of each soil type') |
---|
697 | CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., soiltype, "gather", nbp_glo, index_g) |
---|
698 | ! |
---|
699 | var_name= 'clay_frac' |
---|
700 | CALL ioconf_setatt('UNITS', '-') |
---|
701 | CALL ioconf_setatt('LONG_NAME','Fraction of clay in each mesh') |
---|
702 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g) |
---|
703 | ! |
---|
704 | var_name= 'lai' |
---|
705 | CALL ioconf_setatt('UNITS', '-') |
---|
706 | CALL ioconf_setatt('LONG_NAME','Leaf area index') |
---|
707 | CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g) |
---|
708 | ! |
---|
709 | ! The height of the vegetation could in principle be recalculated at the beginning of the run. |
---|
710 | ! However, this is very tedious, as many special cases have to be taken into account. This variable |
---|
711 | ! is therefore saved in the restart file. |
---|
712 | var_name= 'height' |
---|
713 | CALL ioconf_setatt('UNITS', 'm') |
---|
714 | CALL ioconf_setatt('LONG_NAME','Height of vegetation') |
---|
715 | CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g) |
---|
716 | |
---|
717 | ! |
---|
718 | IF (read_lai)THEN |
---|
719 | ! |
---|
720 | ALLOCATE (laimap(kjpindex,nvm,12),stat=ier) |
---|
721 | IF (ier.NE.0) THEN |
---|
722 | WRITE (numout,*) ' error in laimap allocation. We stop. We need kjpindex*nvm*12 words = ',kjpindex*nvm*12 |
---|
723 | STOP 'slowproc_init' |
---|
724 | END IF |
---|
725 | laimap(:,:,:) = val_exp |
---|
726 | ! |
---|
727 | var_name= 'laimap' |
---|
728 | CALL ioconf_setatt('UNITS', '-') |
---|
729 | CALL ioconf_setatt('LONG_NAME','Leaf area index read') |
---|
730 | CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap) |
---|
731 | ! |
---|
732 | ENDIF |
---|
733 | ! |
---|
734 | !Config Key = SECHIBA_ZCANOP |
---|
735 | !Config Desc = Soil level (m) used for canopy development (if STOMATE disactivated) |
---|
736 | !Config Def = 0.5 |
---|
737 | !Config Help = The temperature at this soil depth is used to determine the LAI when |
---|
738 | !Config STOMATE is not activated. |
---|
739 | ! |
---|
740 | zcanop = 0.5_r_std |
---|
741 | CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std) |
---|
742 | ! |
---|
743 | !Config Key = HYDROL_SOIL_DEPTH |
---|
744 | !Config Desc = Total depth of soil reservoir |
---|
745 | !Config Def = 2. |
---|
746 | ! |
---|
747 | dpu_cste=2. |
---|
748 | CALL getin_p ("HYDROL_SOIL_DEPTH", dpu_cste) |
---|
749 | dpu(:)=dpu_cste |
---|
750 | ! |
---|
751 | !Config Key = HYDROL_HUMCSTE |
---|
752 | !Config Desc = Root profile |
---|
753 | !Config Def = 5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4. |
---|
754 | !Config Help = Default values were defined for 2 meters soil depth. |
---|
755 | !Config For 4 meters soil depth, you may use those ones : |
---|
756 | !Config 5., .4, .4, 1., .8, .8, 1., 1., .8, 4., 1., 4., 1. |
---|
757 | ! |
---|
758 | humcste(:)= & |
---|
759 | & (/5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4./) |
---|
760 | CALL getin_p ("HYDROL_HUMCSTE", humcste) |
---|
761 | |
---|
762 | !MM, T. d'O. : before in constantes_soil : |
---|
763 | ! diaglev = & |
---|
764 | ! & (/ 0.001, 0.004, 0.01, 0.018, 0.045, & |
---|
765 | ! & 0.092, 0.187, 0.375, 0.750, 1.5, & |
---|
766 | ! & 2.0 /) |
---|
767 | ! Place here because it is used in sechiba and stomate (then in teststomate) |
---|
768 | ! to be in sechiba when teststomate will have disapeared. |
---|
769 | !MM Problem here with dpu which depends on soil type |
---|
770 | DO jv = 1, nbdl-1 |
---|
771 | ! first 2.0 is dpu |
---|
772 | ! second 2.0 is average |
---|
773 | diaglev(jv) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / 2.0 |
---|
774 | ENDDO |
---|
775 | diaglev(nbdl) = dpu_cste |
---|
776 | |
---|
777 | ! depth at center of the levels |
---|
778 | zsoil(1) = diaglev(1) / 2. |
---|
779 | DO l = 2, nbdl |
---|
780 | zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2. |
---|
781 | ENDDO |
---|
782 | |
---|
783 | ! index of this level |
---|
784 | vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) ) |
---|
785 | lcanop = vtmp(1) |
---|
786 | |
---|
787 | ! |
---|
788 | ! Interception reservoir coefficient |
---|
789 | ! |
---|
790 | !Config Key = 'SECHIBA_QSINT' |
---|
791 | !Config Desc = Interception reservoir coefficient |
---|
792 | !Config Def = 0.1 |
---|
793 | !Config Help = Transforms leaf area index into size of interception reservoir |
---|
794 | !Config for slowproc_derivvar or stomate |
---|
795 | |
---|
796 | qsintcst = 0.1 |
---|
797 | CALL getin_p('SECHIBA_QSINT', qsintcst) |
---|
798 | WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst |
---|
799 | |
---|
800 | ! |
---|
801 | ! Time step of STOMATE and LAI update |
---|
802 | ! |
---|
803 | !Config Key = DT_SLOW |
---|
804 | !Config Desc = Time step of STOMATE and other slow processes |
---|
805 | !Config Def = 86400. |
---|
806 | !Config Help = Time step (s) of regular update of vegetation |
---|
807 | !Config cover, LAI etc. This is also the time step |
---|
808 | !Config of STOMATE. |
---|
809 | |
---|
810 | dt_slow = one_day |
---|
811 | CALL getin_p('DT_SLOW', dt_slow) |
---|
812 | ! |
---|
813 | |
---|
814 | !Config Key = SLOWPROC_LAI_TEMPDIAG |
---|
815 | !Config Desc = Temperature used for the initial guess of LAI |
---|
816 | !Config Def = 280. |
---|
817 | !Config Help = If there is no LAI in the restart file, we may need |
---|
818 | !Config a temperature that is used to guess the initial LAI. |
---|
819 | ! |
---|
820 | stempdiag_bid = 280. |
---|
821 | CALL getin_p('SLOWPROC_LAI_TEMPDIAG',stempdiag_bid) |
---|
822 | ! |
---|
823 | ! |
---|
824 | ! get restart value if none were found in the restart file |
---|
825 | ! |
---|
826 | !Config Key = AGRICULTURE |
---|
827 | !Config Desc = agriculture allowed? |
---|
828 | !Config Def = y |
---|
829 | !Config Help = With this variable, you can determine |
---|
830 | !Config whether agriculture is allowed |
---|
831 | ! |
---|
832 | agriculture = .TRUE. |
---|
833 | CALL getin_p('AGRICULTURE', agriculture) |
---|
834 | IF ( .NOT. agriculture .AND. land_use ) THEN |
---|
835 | CALL ipslerr (2,'slowproc_init', & |
---|
836 | & 'Problem with agriculture desactivated and Land Use activated.',& |
---|
837 | & 'Are you sure ?', & |
---|
838 | & '(check your parameters).') |
---|
839 | ENDIF |
---|
840 | |
---|
841 | ! |
---|
842 | !Config Key = IMPOSE_VEG |
---|
843 | !Config Desc = Should the vegetation be prescribed |
---|
844 | !Config Def = n |
---|
845 | !Config Help = This flag allows the user to impose a vegetation distribution |
---|
846 | !Config and its characterisitcs. It is espacially interesting for 0D |
---|
847 | !Config simulations. On the globe it does not make too much sense as |
---|
848 | !Config it imposes the same vegetation everywhere |
---|
849 | ! |
---|
850 | impveg = .FALSE. |
---|
851 | CALL getin_p('IMPOSE_VEG', impveg) |
---|
852 | ! |
---|
853 | IF ( impveg ) THEN |
---|
854 | ! |
---|
855 | ! We are on a point and thus we can read the information from the run.def |
---|
856 | ! |
---|
857 | !Config Key = SECHIBA_VEG |
---|
858 | !Config Desc = Vegetation distribution within the mesh (0-dim mode) |
---|
859 | !Config If = IMPOSE_VEG |
---|
860 | !Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0 |
---|
861 | !Config Help = The fraction of vegetation is read from the restart file. If |
---|
862 | !Config it is not found there we will use the values provided here. |
---|
863 | ! |
---|
864 | CALL setvar_p (veget, val_exp, 'SECHIBA_VEG', veget_ori_fixed_test_1) |
---|
865 | |
---|
866 | ! |
---|
867 | !Config Key = SECHIBA_VEGMAX |
---|
868 | !Config Desc = Maximum vegetation distribution within the mesh (0-dim mode) |
---|
869 | !Config If = IMPOSE_VEG |
---|
870 | !Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0 |
---|
871 | !Config Help = The fraction of vegetation is read from the restart file. If |
---|
872 | !Config it is not found there we will use the values provided here. |
---|
873 | ! |
---|
874 | CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1) |
---|
875 | |
---|
876 | ! |
---|
877 | !Config Key = SECHIBA_FRAC_NOBIO |
---|
878 | !Config Desc = Fraction of other surface types within the mesh (0-dim mode) |
---|
879 | !Config If = IMPOSE_VEG |
---|
880 | !Config Def = 0.0 |
---|
881 | !Config Help = The fraction of ice, lakes, etc. is read from the restart file. If |
---|
882 | !Config it is not found there we will use the values provided here. |
---|
883 | !Config For the moment, there is only ice. |
---|
884 | ! |
---|
885 | ! laisser ca tant qu'il n'y a que la glace. Pb avec setvar? |
---|
886 | frac_nobio1 = frac_nobio(1,1) |
---|
887 | CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1) |
---|
888 | frac_nobio(:,:) = frac_nobio1 |
---|
889 | ! CALL setvar (frac_nobio, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1) |
---|
890 | |
---|
891 | ! |
---|
892 | !Config Key = SECHIBA_LAI |
---|
893 | !Config Desc = LAI for all vegetation types (0-dim mode) |
---|
894 | !Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2. |
---|
895 | !Config If = IMPOSE_VEG |
---|
896 | !Config Help = The maximum LAI used in the 0dim mode. The values should be found |
---|
897 | !Config in the restart file. The new values of LAI will be computed anyway |
---|
898 | !Config at the end of the current day. The need for this variable is caused |
---|
899 | !Config by the fact that the model may stop during a day and thus we have not |
---|
900 | !Config yet been through the routines which compute the new surface conditions. |
---|
901 | ! |
---|
902 | CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax) |
---|
903 | |
---|
904 | !Config Key = SOIL_FRACTIONS |
---|
905 | !Config Desc = Fraction of the 3 soil types (0-dim mode) |
---|
906 | !Config Def = 0.28, 0.52, 0.20 |
---|
907 | !Config If = IMPOSE_VEG |
---|
908 | !Config Help = Determines the fraction for the 3 soil types |
---|
909 | !Config in the mesh in the following order : sand loam and clay. |
---|
910 | ! |
---|
911 | CALL setvar_p (soiltype, val_exp, 'SOIL_FRACTIONS', soiltype_default) |
---|
912 | |
---|
913 | !Config Key = CLAY_FRACTION |
---|
914 | !Config Desc = Fraction of the clay fraction (0-dim mode) |
---|
915 | !Config Def = 0.2 |
---|
916 | !Config If = IMPOSE_VEG |
---|
917 | !Config Help = Determines the fraction of clay in the grid box. |
---|
918 | ! |
---|
919 | CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default) |
---|
920 | |
---|
921 | ! |
---|
922 | !Config Key = SLOWPROC_HEIGHT |
---|
923 | !Config Desc = Height for all vegetation types (m) |
---|
924 | !Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0 |
---|
925 | !Config If = IMPOSE_VEG |
---|
926 | !Config Help = The height used in the 0dim mode. The values should be found |
---|
927 | !Config in the restart file. The new values of height will be computed anyway |
---|
928 | !Config at the end of the current day. The need for this variable is caused |
---|
929 | !Config by the fact that the model may stop during a day and thus we have not |
---|
930 | !Config yet been through the routines which compute the new surface conditions. |
---|
931 | ! |
---|
932 | CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc) |
---|
933 | |
---|
934 | ELSE |
---|
935 | ! |
---|
936 | ! We are in the full 2-D case thus we need to do the interpolation if the potential vegetation |
---|
937 | ! is not available |
---|
938 | ! |
---|
939 | IF ( ALL( veget_max(:,:) .EQ. val_exp ) .OR. ALL( frac_nobio(:,:) .EQ. val_exp ) ) THEN |
---|
940 | |
---|
941 | IF ( .NOT. land_use ) THEN |
---|
942 | |
---|
943 | !Config Key = SLOWPROC_VEGET_OLD_INTERPOL |
---|
944 | !Config Desc = Flag to use old "interpolation" of vegetation map. |
---|
945 | !Config If = NOT IMPOSE_VEG and NOT LAND_USE |
---|
946 | !Config Def = FALSE |
---|
947 | !Config Help = If you want to recover the old (ie orchidee_1_2 branch) |
---|
948 | !Config "interpolation" of vegetation map. |
---|
949 | ! |
---|
950 | old_veget = .FALSE. |
---|
951 | CALL getin_p('SLOWPROC_VEGET_OLD_INTERPOL',old_veget) |
---|
952 | |
---|
953 | ! The interpolation of vegetation has changed. |
---|
954 | IF (is_root_prc) THEN |
---|
955 | IF ( .NOT. old_veget ) THEN |
---|
956 | ! NEW slowproc interpol : |
---|
957 | CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, veget_max_g, frac_nobio_g) |
---|
958 | ELSE |
---|
959 | ! OLD slowproc interpol : |
---|
960 | CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, veget_max_g, frac_nobio_g) |
---|
961 | ENDIF |
---|
962 | ENDIF |
---|
963 | CALL scatter(veget_max_g,veget_max) |
---|
964 | CALL scatter(frac_nobio_g, frac_nobio) |
---|
965 | ! |
---|
966 | IF ( control%ok_dgvm ) THEN |
---|
967 | ! |
---|
968 | ! If we are dealing with dynamic vegetation then all |
---|
969 | ! natural PFTs should be set to veget_max = 0 |
---|
970 | ! And in case no agriculture is desired, agriculture PFTS should be |
---|
971 | ! set to 0 as well |
---|
972 | |
---|
973 | IF (agriculture) THEN |
---|
974 | veget_max(:,2:nvm-2) = zero |
---|
975 | DO ji = 1, kjpindex |
---|
976 | veget_max(ji,1) = un - SUM(veget_max(ji,nvm-1:nvm)) & |
---|
977 | - SUM(frac_nobio(ji,:)) |
---|
978 | ENDDO |
---|
979 | ELSE |
---|
980 | veget_max(:,:) = zero |
---|
981 | DO ji = 1, kjpindex |
---|
982 | veget_max(ji,1) = un - SUM(frac_nobio(ji,:)) |
---|
983 | ENDDO |
---|
984 | ENDIF |
---|
985 | ! |
---|
986 | ENDIF |
---|
987 | ELSE |
---|
988 | WRITE(numout,*) 'We initialize land use veget for year =' , veget_year |
---|
989 | ! If restart doesn't contain veget, then it is the first computation |
---|
990 | CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & |
---|
991 | & veget_max, frac_nobio, veget_max, frac_nobio, veget_year, init=.TRUE.) |
---|
992 | ! |
---|
993 | IF ( control%ok_dgvm ) THEN |
---|
994 | ! |
---|
995 | ! If we are dealing with dynamic vegetation then all |
---|
996 | ! natural PFTs should be set to veget_max = 0 |
---|
997 | ! And in case no agriculture is desired, agriculture PFTS should be |
---|
998 | ! set to 0 as well |
---|
999 | |
---|
1000 | IF (agriculture) THEN |
---|
1001 | veget_max(:,2:nvm-2) = zero |
---|
1002 | DO ji = 1, kjpindex |
---|
1003 | veget_max(ji,1) = un - SUM(veget_max(ji,nvm-1:nvm)) & |
---|
1004 | - SUM(frac_nobio(ji,:)) |
---|
1005 | ENDDO |
---|
1006 | ELSE |
---|
1007 | veget_max(:,:) = zero |
---|
1008 | DO ji = 1, kjpindex |
---|
1009 | veget_max(ji,1) = un - SUM(frac_nobio(ji,:)) |
---|
1010 | ENDDO |
---|
1011 | ENDIF |
---|
1012 | ! |
---|
1013 | ENDIF |
---|
1014 | ! |
---|
1015 | ENDIF |
---|
1016 | ! |
---|
1017 | ELSE |
---|
1018 | ! WITH restarts for vegetation and DGVM and NO AGRICULTURE |
---|
1019 | IF ( control%ok_dgvm ) THEN |
---|
1020 | ! |
---|
1021 | ! If we are dealing with dynamic vegetation then all |
---|
1022 | ! natural PFTs should be set to veget_max = 0 |
---|
1023 | ! And in case no agriculture is desired, agriculture PFTS should be |
---|
1024 | ! set to 0 as well |
---|
1025 | ! |
---|
1026 | IF (.NOT. agriculture) THEN |
---|
1027 | DO ji = 1, kjpindex |
---|
1028 | veget_max(ji,1) = veget_max(ji,1) + SUM(veget_max(ji,nvm-1:nvm)) |
---|
1029 | ENDDO |
---|
1030 | veget_max(ji,nvm-1:nvm) = zero |
---|
1031 | ENDIF |
---|
1032 | ! |
---|
1033 | ENDIF |
---|
1034 | ! |
---|
1035 | ENDIF |
---|
1036 | ! |
---|
1037 | totfrac_nobio(:) = zero |
---|
1038 | DO j = 1, nnobio |
---|
1039 | totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,j) |
---|
1040 | ENDDO |
---|
1041 | ! |
---|
1042 | ! |
---|
1043 | IF (read_lai) THEN |
---|
1044 | |
---|
1045 | !Config Key = SLOWPROC_LAI_OLD_INTERPOL |
---|
1046 | !Config Desc = Flag to use old "interpolation" of LAI |
---|
1047 | !Config If = LAI_MAP |
---|
1048 | !Config Def = FALSE |
---|
1049 | !Config Help = If you want to recover the old (ie orchidee_1_2 branch) |
---|
1050 | !Config "interpolation" of LAI map. |
---|
1051 | ! |
---|
1052 | old_lai = .FALSE. |
---|
1053 | CALL getin_p('SLOWPROC_LAI_OLD_INTERPOL',old_lai) |
---|
1054 | |
---|
1055 | ! |
---|
1056 | ! In case the LAI map was not found in the restart then we need to recompute it |
---|
1057 | ! |
---|
1058 | IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN |
---|
1059 | ! The interpolation of LAI has changed. |
---|
1060 | IF ( .NOT. old_lai ) THEN |
---|
1061 | ! NEW slowproc interlai : |
---|
1062 | CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap) |
---|
1063 | ELSE |
---|
1064 | ! OLD slowproc interlai : |
---|
1065 | CALL slowproc_interlai(kjpindex, lalo, resolution, laimap) |
---|
1066 | ENDIF |
---|
1067 | ! |
---|
1068 | ! Compute the current LAI just to be sure |
---|
1069 | ! |
---|
1070 | stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid |
---|
1071 | CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, & |
---|
1072 | lalo,resolution,lai,month,day,read_lai,laimap) |
---|
1073 | ! |
---|
1074 | ! Make sure that we redo the computation when we will be back in slowproc_main |
---|
1075 | day_counter = dt_slow - dtradia |
---|
1076 | ! |
---|
1077 | ENDIF |
---|
1078 | ! |
---|
1079 | ENDIF |
---|
1080 | ! |
---|
1081 | IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN |
---|
1082 | ! |
---|
1083 | ! Get a first guess at the LAI |
---|
1084 | ! |
---|
1085 | IF ( read_lai ) THEN |
---|
1086 | IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN |
---|
1087 | ! The interpolation of LAI has changed. |
---|
1088 | IF ( .NOT. old_lai ) THEN |
---|
1089 | ! NEW slowproc interlai : |
---|
1090 | CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap) |
---|
1091 | ELSE |
---|
1092 | ! OLD slowproc interlai : |
---|
1093 | CALL slowproc_interlai(kjpindex, lalo, resolution, laimap) |
---|
1094 | ENDIF |
---|
1095 | ENDIF |
---|
1096 | ENDIF |
---|
1097 | ! |
---|
1098 | stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid |
---|
1099 | CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, & |
---|
1100 | lalo,resolution,lai,month,day,read_lai,laimap) |
---|
1101 | |
---|
1102 | ! Make sure that we redo the computation when we will be back in slowproc_main |
---|
1103 | day_counter = dt_slow - dtradia |
---|
1104 | |
---|
1105 | ENDIF |
---|
1106 | |
---|
1107 | IF ( MINVAL(veget) .EQ. MAXVAL(veget) .AND. MAXVAL(veget) .EQ. val_exp) THEN |
---|
1108 | |
---|
1109 | ! Impose height |
---|
1110 | DO jv = 1, nvm |
---|
1111 | height(:,jv) = height_presc(jv) |
---|
1112 | ENDDO |
---|
1113 | |
---|
1114 | ! Have a first guess at the vegetation fraction |
---|
1115 | CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) |
---|
1116 | |
---|
1117 | ENDIF |
---|
1118 | |
---|
1119 | IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.& |
---|
1120 | & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN |
---|
1121 | |
---|
1122 | CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction) |
---|
1123 | |
---|
1124 | ENDIF |
---|
1125 | |
---|
1126 | ENDIF |
---|
1127 | ! |
---|
1128 | ! Select the preferences for the distribution of the 13 PFTs on the 3 soil types. |
---|
1129 | ! |
---|
1130 | vegsoil_dist='EQUI' |
---|
1131 | ! |
---|
1132 | SELECT CASE(vegsoil_dist) |
---|
1133 | ! |
---|
1134 | ! A reasonable choice |
---|
1135 | ! |
---|
1136 | CASE('MAXR') |
---|
1137 | pref_soil_veg(:,1) = (/ 1, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 /) |
---|
1138 | pref_soil_veg(:,2) = (/ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 /) |
---|
1139 | pref_soil_veg(:,3) = (/ 3, 1, 1, 1, 1, 1 ,1 ,1 ,1 ,1 ,1 ,1, 1 /) |
---|
1140 | ! |
---|
1141 | ! Current default : equidistribution. |
---|
1142 | ! |
---|
1143 | CASE('EQUI') |
---|
1144 | ! |
---|
1145 | pref_soil_veg(:,:) = zero |
---|
1146 | ! |
---|
1147 | CASE DEFAULT |
---|
1148 | ! |
---|
1149 | WRITE(numout,*) 'The vegetation soil type preference you have chosen does not exist' |
---|
1150 | WRITE(numout,*) 'You chose :', vegsoil_dist |
---|
1151 | STOP 'slowproc_init' |
---|
1152 | ! |
---|
1153 | END SELECT |
---|
1154 | ! |
---|
1155 | ! |
---|
1156 | ! calculate total fraction of the mesh that is covered by particular surface types: ice, lakes, ... |
---|
1157 | ! |
---|
1158 | totfrac_nobio(:) = zero |
---|
1159 | DO jv = 1, nnobio |
---|
1160 | totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv) |
---|
1161 | ENDDO |
---|
1162 | |
---|
1163 | l_first_slowproc = .FALSE. |
---|
1164 | |
---|
1165 | IF (long_print) WRITE (numout,*) ' slowproc_init done ' |
---|
1166 | |
---|
1167 | END SUBROUTINE slowproc_init |
---|
1168 | !! |
---|
1169 | !! Clear Memory |
---|
1170 | !! |
---|
1171 | SUBROUTINE slowproc_clear |
---|
1172 | |
---|
1173 | |
---|
1174 | l_first_slowproc = .TRUE. |
---|
1175 | IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction) |
---|
1176 | IF (ALLOCATED (laimap)) DEALLOCATE (laimap) |
---|
1177 | ! |
---|
1178 | ! |
---|
1179 | IF (ALLOCATED (veget_nextyear)) DEALLOCATE (veget_nextyear) |
---|
1180 | IF (ALLOCATED (frac_nobio_nextyear)) DEALLOCATE (frac_nobio_nextyear) |
---|
1181 | IF (ALLOCATED (totfrac_nobio_nextyear)) DEALLOCATE (totfrac_nobio_nextyear) |
---|
1182 | ! |
---|
1183 | CALL stomate_clear |
---|
1184 | ! |
---|
1185 | END SUBROUTINE slowproc_clear |
---|
1186 | !! |
---|
1187 | !! Derive some variables |
---|
1188 | !! |
---|
1189 | SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, & |
---|
1190 | qsintmax, deadleaf_cover, assim_param, height) |
---|
1191 | |
---|
1192 | ! interface description |
---|
1193 | ! input scalar |
---|
1194 | INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size |
---|
1195 | ! input fields |
---|
1196 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type |
---|
1197 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Surface foliere |
---|
1198 | ! output scalar |
---|
1199 | ! output fields |
---|
1200 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception |
---|
1201 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover!! fraction of soil covered by dead leaves |
---|
1202 | REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temps & vmax for photosynthesis |
---|
1203 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height (m) |
---|
1204 | ! |
---|
1205 | ! local declaration |
---|
1206 | ! |
---|
1207 | INTEGER(i_std) :: ji, jv |
---|
1208 | |
---|
1209 | ! |
---|
1210 | ! 1 Assimilation parameters |
---|
1211 | ! |
---|
1212 | DO jv = 1, nvm |
---|
1213 | assim_param(:,jv,itmin) = co2_tmin_fix(jv) + tp_00 |
---|
1214 | assim_param(:,jv,itopt) = co2_topt_fix(jv) + tp_00 |
---|
1215 | assim_param(:,jv,itmax) = co2_tmax_fix(jv) + tp_00 |
---|
1216 | assim_param(:,jv,ivcmax) = vcmax_fix(jv) |
---|
1217 | assim_param(:,jv,ivjmax) = vjmax_fix(jv) |
---|
1218 | ENDDO |
---|
1219 | |
---|
1220 | ! |
---|
1221 | ! 2 fraction of soil covered by dead leaves |
---|
1222 | ! |
---|
1223 | deadleaf_cover(:) = zero |
---|
1224 | |
---|
1225 | ! |
---|
1226 | ! 3 height |
---|
1227 | ! |
---|
1228 | DO jv = 1, nvm |
---|
1229 | height(:,jv) = height_presc(jv) |
---|
1230 | ENDDO |
---|
1231 | ! |
---|
1232 | ! 4 qsintmax |
---|
1233 | ! |
---|
1234 | qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:) |
---|
1235 | ! Ajout Nathalie - Juillet 2006 |
---|
1236 | qsintmax(:,1) = zero |
---|
1237 | |
---|
1238 | END SUBROUTINE slowproc_derivvar |
---|
1239 | !! |
---|
1240 | !! |
---|
1241 | !! |
---|
1242 | SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean) |
---|
1243 | |
---|
1244 | ! Accumulates field_in over a period of dt_tot. |
---|
1245 | ! Has to be called at every time step (dt). |
---|
1246 | ! Mean value is calculated if ldmean=.TRUE. |
---|
1247 | ! field_mean must be initialized outside of this routine! |
---|
1248 | |
---|
1249 | ! |
---|
1250 | ! 0 declarations |
---|
1251 | ! |
---|
1252 | |
---|
1253 | ! 0.1 input |
---|
1254 | |
---|
1255 | ! Domain size |
---|
1256 | INTEGER(i_std), INTENT(in) :: npts |
---|
1257 | ! 2nd dimension (1 or npft) |
---|
1258 | INTEGER(i_std), INTENT(in) :: n_dim2 |
---|
1259 | ! Time step of STOMATE (days) |
---|
1260 | REAL(r_std), INTENT(in) :: dt_tot |
---|
1261 | ! Time step in days |
---|
1262 | REAL(r_std), INTENT(in) :: dt |
---|
1263 | ! Calculate mean ? |
---|
1264 | LOGICAL, INTENT(in) :: ldmean |
---|
1265 | ! Daily field |
---|
1266 | REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in |
---|
1267 | |
---|
1268 | ! 0.2 modified field |
---|
1269 | |
---|
1270 | ! Annual field |
---|
1271 | REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean |
---|
1272 | |
---|
1273 | ! ========================================================================= |
---|
1274 | |
---|
1275 | ! |
---|
1276 | ! 1 accumulation |
---|
1277 | ! |
---|
1278 | |
---|
1279 | field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt |
---|
1280 | |
---|
1281 | ! |
---|
1282 | ! 2 mean fields |
---|
1283 | ! |
---|
1284 | |
---|
1285 | IF (ldmean) THEN |
---|
1286 | field_mean(:,:) = field_mean(:,:) / dt_tot |
---|
1287 | ENDIF |
---|
1288 | |
---|
1289 | END SUBROUTINE slowproc_mean |
---|
1290 | !! |
---|
1291 | !! |
---|
1292 | !! |
---|
1293 | SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long) |
---|
1294 | |
---|
1295 | ! Calculates a temporally smoothed field (field_long) from instantaneous |
---|
1296 | ! input fields. |
---|
1297 | ! Time constant tau determines the strength of the smoothing. |
---|
1298 | ! For tau -> infty, field_long becomes the true mean value of field_inst (but |
---|
1299 | ! the spinup becomes infinietly long, too). |
---|
1300 | ! field_long must be initialized outside of this routine! |
---|
1301 | |
---|
1302 | ! |
---|
1303 | ! 0 declarations |
---|
1304 | ! |
---|
1305 | |
---|
1306 | ! 0.1 input |
---|
1307 | |
---|
1308 | ! Domain size |
---|
1309 | INTEGER(i_std), INTENT(in) :: npts |
---|
1310 | ! 2nd dimension (1 or npft) |
---|
1311 | INTEGER(i_std), INTENT(in) :: n_dim2 |
---|
1312 | ! Time step |
---|
1313 | REAL(r_std), INTENT(in) :: dt |
---|
1314 | ! Integration time constant (has to have same unit as dt!) |
---|
1315 | REAL(r_std), INTENT(in) :: tau |
---|
1316 | ! Instantaneous field |
---|
1317 | REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst |
---|
1318 | |
---|
1319 | ! 0.2 modified field |
---|
1320 | |
---|
1321 | ! Long-term field |
---|
1322 | REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long |
---|
1323 | |
---|
1324 | ! ========================================================================= |
---|
1325 | |
---|
1326 | ! |
---|
1327 | ! 1 test coherence |
---|
1328 | ! |
---|
1329 | |
---|
1330 | IF ( ( tau .LT. dt ) .OR. ( dt .LE. 0. ) .OR. ( tau .LE. 0. ) ) THEN |
---|
1331 | WRITE(numout,*) 'slowproc_long: Problem with time steps' |
---|
1332 | WRITE(numout,*) 'dt=',dt |
---|
1333 | WRITE(numout,*) 'tau=',tau |
---|
1334 | ENDIF |
---|
1335 | |
---|
1336 | ! |
---|
1337 | ! 2 integration |
---|
1338 | ! |
---|
1339 | |
---|
1340 | field_long(:,:) = ( field_inst(:,:)*dt + field_long(:,:)*(tau-dt) ) / tau |
---|
1341 | |
---|
1342 | END SUBROUTINE slowproc_long |
---|
1343 | !! |
---|
1344 | !! |
---|
1345 | !! |
---|
1346 | SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) |
---|
1347 | ! |
---|
1348 | ! 0. Declarations |
---|
1349 | ! |
---|
1350 | |
---|
1351 | ! 0.1 Input |
---|
1352 | INTEGER(i_std), INTENT(in) :: kjpindex |
---|
1353 | REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai |
---|
1354 | |
---|
1355 | ! 0.2 Modified |
---|
1356 | REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio |
---|
1357 | REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max |
---|
1358 | |
---|
1359 | ! 0.3 Output |
---|
1360 | REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget |
---|
1361 | |
---|
1362 | ! 0.4 Local |
---|
1363 | REAL(r_std), DIMENSION(kjpindex) :: fracsum |
---|
1364 | INTEGER(i_std) :: nbad |
---|
1365 | INTEGER(i_std) :: ji, jv |
---|
1366 | ! Test Nathalie |
---|
1367 | REAL(r_std) :: SUMveg |
---|
1368 | !!$ REAL(r_std) :: trans_veg |
---|
1369 | |
---|
1370 | ! |
---|
1371 | ! 1. Sum up veget_max and frac_nobio and test if sum is equal to 1 |
---|
1372 | ! |
---|
1373 | ! |
---|
1374 | ! 1.1 Sum up |
---|
1375 | ! |
---|
1376 | fracsum(:) = 0. |
---|
1377 | DO jv = 1, nnobio |
---|
1378 | DO ji = 1, kjpindex |
---|
1379 | fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv) |
---|
1380 | ENDDO |
---|
1381 | ENDDO |
---|
1382 | DO jv = 1, nvm |
---|
1383 | DO ji = 1, kjpindex |
---|
1384 | fracsum(ji) = fracsum(ji) + veget_max(ji,jv) |
---|
1385 | ENDDO |
---|
1386 | ENDDO |
---|
1387 | ! |
---|
1388 | ! 1.2 stop if there is a severe problem, that is an error above 0.01% |
---|
1389 | ! The rest will be scaled |
---|
1390 | ! |
---|
1391 | nbad = COUNT( ABS(fracsum(:)-un) .GT. 0.0001 ) |
---|
1392 | IF ( nbad .GT. 0 ) THEN |
---|
1393 | WRITE(numout,*) 'Problem with total surface areas.' |
---|
1394 | DO ji = 1, kjpindex |
---|
1395 | IF ( ABS(fracsum(ji)-un) .GT. 0.0001 ) THEN |
---|
1396 | WRITE(numout,*) 'Point :', ji |
---|
1397 | WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:) |
---|
1398 | WRITE(numout,*) ' Veget_max :', veget_max(ji,:) |
---|
1399 | WRITE(numout,*) ' Fracsum :', fracsum(ji), EPSILON(un) |
---|
1400 | WRITE(numout,*) ' The error is :', un - ( SUM(frac_nobio(ji,:)) + SUM(veget_max(ji,:)) ) |
---|
1401 | STOP 'slowproc_veget' |
---|
1402 | ENDIF |
---|
1403 | ENDDO |
---|
1404 | ENDIF |
---|
1405 | ! |
---|
1406 | ! 1.3 correction at places where the problem is surely precision-related |
---|
1407 | ! |
---|
1408 | nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) ) |
---|
1409 | ! |
---|
1410 | IF ( nbad .GT. 0 ) THEN |
---|
1411 | ! |
---|
1412 | IF ( long_print ) THEN |
---|
1413 | WRITE(numout,*) 'WARNING! scaling frac_nobio and veget_max at', nbad, ' points' |
---|
1414 | ENDIF |
---|
1415 | ! |
---|
1416 | DO jv = 1, nnobio |
---|
1417 | DO ji = 1, kjpindex |
---|
1418 | IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN |
---|
1419 | frac_nobio(ji,jv) = frac_nobio(ji,jv)/fracsum(ji) |
---|
1420 | ENDIF |
---|
1421 | ENDDO |
---|
1422 | ENDDO |
---|
1423 | ! |
---|
1424 | DO jv = 1, nvm |
---|
1425 | DO ji = 1, kjpindex |
---|
1426 | IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN |
---|
1427 | veget_max(ji,jv) = veget_max(ji,jv)/fracsum(ji) |
---|
1428 | ENDIF |
---|
1429 | ENDDO |
---|
1430 | ENDDO |
---|
1431 | ! |
---|
1432 | ENDIF |
---|
1433 | |
---|
1434 | ! |
---|
1435 | ! 2. Set veget equal to veget_max |
---|
1436 | ! |
---|
1437 | DO jv = 1, nvm |
---|
1438 | DO ji = 1, kjpindex |
---|
1439 | veget(ji,jv) = veget_max(ji,jv) |
---|
1440 | ENDDO |
---|
1441 | ENDDO |
---|
1442 | ! |
---|
1443 | ! 3. if lai of a vegetation type (jv > 1) is small, increase soil part |
---|
1444 | ! |
---|
1445 | ! Nathalie - nouveau calcul de veget |
---|
1446 | !!$ DO jv = 1,nvm |
---|
1447 | !!$ DO ji = 1, kjpindex |
---|
1448 | !!$ |
---|
1449 | !!$ IF ((jv .GT. 1) .AND. (lai(ji,jv) .LT. 0.5)) THEN |
---|
1450 | !!$ trans_veg = 2.*(0.5 - lai(ji,jv)) * veget(ji,jv) |
---|
1451 | !!$ ! We limit the smallest fraction to 0.5% |
---|
1452 | !!$ IF ( veget(ji,jv) - trans_veg .GT. min_vegfrac ) THEN |
---|
1453 | !!$ veget(ji,1) = veget(ji,1) + trans_veg |
---|
1454 | !!$ veget(ji,jv) = veget(ji,jv) - trans_veg |
---|
1455 | !!$ ELSE |
---|
1456 | !!$ veget(ji,1) = veget(ji,1) + veget(ji,jv) |
---|
1457 | !!$ veget(ji,jv) = zero |
---|
1458 | !!$ ENDIF |
---|
1459 | !!$ ENDIF |
---|
1460 | !!$ |
---|
1461 | !!$ ENDDO |
---|
1462 | !!$ ENDDO |
---|
1463 | ! Ajout Nouveau calcul (stomate-like) |
---|
1464 | DO ji = 1, kjpindex |
---|
1465 | SUMveg = 0.0 |
---|
1466 | veget(ji,1) = veget_max(ji,1) |
---|
1467 | DO jv = 2, nvm |
---|
1468 | veget(ji,jv) = veget_max(ji,jv) * ( 1. - exp( - lai(ji,jv) * ext_coef(jv) ) ) |
---|
1469 | veget(ji,1) = veget(ji,1) + (veget_max(ji,jv) - veget(ji,jv)) |
---|
1470 | SUMveg = SUMveg + veget(ji,jv) |
---|
1471 | ENDDO |
---|
1472 | SUMveg = SUMveg + veget(ji,1) + SUM(frac_nobio(ji,:)) |
---|
1473 | IF (SUMveg .LT. 0.99999) THEN |
---|
1474 | WRITE(numout,*)' ATTENTION, en ji, SUMveg LT 1: ', ji, SUMveg |
---|
1475 | WRITE(numout,*)' frac_nobio = ',SUM(frac_nobio(ji,:)) |
---|
1476 | WRITE(numout,*) ' ',veget(ji,:) |
---|
1477 | ENDIF |
---|
1478 | ENDDO |
---|
1479 | ! |
---|
1480 | ! 4. Sum up surface fractions and test if the sum is equal to 1 |
---|
1481 | ! |
---|
1482 | |
---|
1483 | ! |
---|
1484 | ! 4.1 Sum up |
---|
1485 | ! |
---|
1486 | fracsum(:) = 0. |
---|
1487 | DO jv = 1, nnobio |
---|
1488 | DO ji = 1, kjpindex |
---|
1489 | fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv) |
---|
1490 | ENDDO |
---|
1491 | ENDDO |
---|
1492 | DO jv = 1, nvm |
---|
1493 | DO ji = 1, kjpindex |
---|
1494 | fracsum(ji) = fracsum(ji) + veget_max(ji,jv) |
---|
1495 | ENDDO |
---|
1496 | ENDDO |
---|
1497 | |
---|
1498 | ! |
---|
1499 | ! 4.2 stop if there is a severe problem |
---|
1500 | ! |
---|
1501 | nbad = COUNT( ABS(fracsum(:)-un) .GT. (REAL(nvm+nnobio,r_std)*EPSILON(un)) ) |
---|
1502 | ! |
---|
1503 | IF ( nbad .GT. 0 ) THEN |
---|
1504 | WRITE(numout,*) 'Problem with veget or frac_nobio.' |
---|
1505 | DO ji = 1, kjpindex |
---|
1506 | IF ( ABS(fracsum(ji)-un) .GT. (10.*EPSILON(un)) ) THEN |
---|
1507 | WRITE(numout,*) 'Point :', ji |
---|
1508 | WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:) |
---|
1509 | WRITE(numout,*) ' Veget :', veget(ji,:) |
---|
1510 | WRITE(numout,*) ' The error is :', un - (SUM(frac_nobio(ji,:)) + SUM(veget(ji,:))) |
---|
1511 | STOP 'slowproc_veget' |
---|
1512 | ENDIF |
---|
1513 | ENDDO |
---|
1514 | ENDIF |
---|
1515 | |
---|
1516 | ! |
---|
1517 | ! 4.3 correction at places where the problem is surely precision-related |
---|
1518 | ! |
---|
1519 | nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) ) |
---|
1520 | ! |
---|
1521 | IF ( nbad .GT. 0 ) THEN |
---|
1522 | ! |
---|
1523 | IF ( long_print ) THEN |
---|
1524 | WRITE(numout,*) 'slowproc_veget: scaling veget at', nbad, ' points' |
---|
1525 | ENDIF |
---|
1526 | ! |
---|
1527 | DO jv = 1, nvm |
---|
1528 | DO ji = 1, kjpindex |
---|
1529 | IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN |
---|
1530 | veget(ji,jv) = veget(ji,jv) / fracsum(ji) |
---|
1531 | ENDIF |
---|
1532 | ENDDO |
---|
1533 | ENDDO |
---|
1534 | ! |
---|
1535 | DO jv = 1, nnobio |
---|
1536 | DO ji = 1, kjpindex |
---|
1537 | IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN |
---|
1538 | frac_nobio(ji,jv) = frac_nobio(ji,jv) / fracsum(ji) |
---|
1539 | ENDIF |
---|
1540 | ENDDO |
---|
1541 | ENDDO |
---|
1542 | ! |
---|
1543 | ENDIF |
---|
1544 | |
---|
1545 | END SUBROUTINE slowproc_veget |
---|
1546 | !! |
---|
1547 | !! |
---|
1548 | !! |
---|
1549 | SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,mm,dd,read_lai,laimap) |
---|
1550 | ! |
---|
1551 | ! 0. Declarations |
---|
1552 | ! |
---|
1553 | |
---|
1554 | ! 0.1 Input |
---|
1555 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size |
---|
1556 | INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI |
---|
1557 | INTEGER(i_std), INTENT(in) :: mm, dd |
---|
1558 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature |
---|
1559 | REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) |
---|
1560 | REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) |
---|
1561 | |
---|
1562 | REAL(r_std), DIMENSION(kjpindex,nvm,12), INTENT(in) :: laimap !! LAI lue |
---|
1563 | LOGICAL, INTENT(in) :: read_lai |
---|
1564 | ! 0.2 Output |
---|
1565 | REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! LAI |
---|
1566 | |
---|
1567 | ! 0.3 Local |
---|
1568 | INTEGER(i_std) :: ji,jv |
---|
1569 | |
---|
1570 | |
---|
1571 | ! Test Nathalie. On impose LAI PFT 1 a 0 |
---|
1572 | ! On boucle sur 2,nvm au lieu de 1,nvm |
---|
1573 | lai(: ,1) = 0.0 |
---|
1574 | DO jv = 2,nvm |
---|
1575 | !!$ DO jv = 1,nvm |
---|
1576 | |
---|
1577 | SELECT CASE (type_of_lai(jv)) |
---|
1578 | |
---|
1579 | CASE ("mean ") |
---|
1580 | ! |
---|
1581 | ! 1. do the interpolation between laimax and laimin |
---|
1582 | ! |
---|
1583 | IF ( .NOT. read_lai ) THEN |
---|
1584 | lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv)) |
---|
1585 | ELSE |
---|
1586 | DO ji = 1,kjpindex |
---|
1587 | lai(ji,jv) = MAXVAL(laimap(ji,jv,:)) |
---|
1588 | ENDDO |
---|
1589 | ENDIF |
---|
1590 | ! |
---|
1591 | CASE ("inter") |
---|
1592 | ! |
---|
1593 | ! 2. do the interpolation between laimax and laimin |
---|
1594 | ! |
---|
1595 | IF ( .NOT. read_lai ) THEN |
---|
1596 | DO ji = 1,kjpindex |
---|
1597 | lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv)) |
---|
1598 | ENDDO |
---|
1599 | ELSE |
---|
1600 | IF (mm .EQ. 1 ) THEN |
---|
1601 | IF (dd .LE. 15) THEN |
---|
1602 | lai(:,jv) = laimap(:,jv,12)*(1-(dd+15)/30.) + laimap(:,jv,1)*((dd+15)/30.) |
---|
1603 | ELSE |
---|
1604 | lai(:,jv) = laimap(:,jv,1)*(1-(dd-15)/30.) + laimap(:,jv,2)*((dd-15)/30.) |
---|
1605 | ENDIF |
---|
1606 | ELSE IF (mm .EQ. 12) THEN |
---|
1607 | IF (dd .LE. 15) THEN |
---|
1608 | lai(:,jv) = laimap(:,jv,11)*(1-(dd+15)/30.) + laimap(:,jv,12)*((dd+15)/30.) |
---|
1609 | ELSE |
---|
1610 | lai(:,jv) = laimap(:,jv,12)*(1-(dd-15)/30.) + laimap(:,jv,1)*((dd-15)/30.) |
---|
1611 | ENDIF |
---|
1612 | ELSE |
---|
1613 | IF (dd .LE. 15) THEN |
---|
1614 | lai(:,jv) = laimap(:,jv,mm-1)*(1-(dd+15)/30.) + laimap(:,jv,mm)*((dd+15)/30.) |
---|
1615 | ELSE |
---|
1616 | lai(:,jv) = laimap(:,jv,mm)*(1-(dd-15)/30.) + laimap(:,jv,mm+1)*((dd-15)/30.) |
---|
1617 | ENDIF |
---|
1618 | ENDIF |
---|
1619 | ENDIF |
---|
1620 | ! |
---|
1621 | CASE default |
---|
1622 | ! |
---|
1623 | ! 3. Problem |
---|
1624 | ! |
---|
1625 | WRITE (numout,*) 'This kind of lai choice is not possible. '// & |
---|
1626 | ' We stop with type_of_lai ',jv,' = ', type_of_lai(jv) |
---|
1627 | STOP 'slowproc_lai' |
---|
1628 | |
---|
1629 | END SELECT |
---|
1630 | |
---|
1631 | ENDDO |
---|
1632 | |
---|
1633 | END SUBROUTINE slowproc_lai |
---|
1634 | !! |
---|
1635 | !! Interpolate the LAI map to the grid of the model |
---|
1636 | !MM TAG 1.6 model ! |
---|
1637 | !! |
---|
1638 | SUBROUTINE slowproc_interlai_OLD(nbpt, lalo, resolution, laimap) |
---|
1639 | ! |
---|
1640 | ! |
---|
1641 | ! |
---|
1642 | ! 0.1 INPUT |
---|
1643 | ! |
---|
1644 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
1645 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
1646 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
1647 | ! |
---|
1648 | ! 0.2 OUTPUT |
---|
1649 | ! |
---|
1650 | REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned |
---|
1651 | ! |
---|
1652 | ! 0.3 LOCAL |
---|
1653 | ! |
---|
1654 | REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8 |
---|
1655 | ! |
---|
1656 | ! |
---|
1657 | CHARACTER(LEN=80) :: filename |
---|
1658 | INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj |
---|
1659 | REAL(r_std) :: lev(1), date, dt, coslat, pi |
---|
1660 | INTEGER(i_std) :: itau(1) |
---|
1661 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu |
---|
1662 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask |
---|
1663 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful |
---|
1664 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: laimaporig |
---|
1665 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu |
---|
1666 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low |
---|
1667 | INTEGER, DIMENSION(nbpt) :: n_origlai |
---|
1668 | INTEGER, DIMENSION(nbpt) :: n_found |
---|
1669 | REAL(r_std), DIMENSION(nbpt,nvm) :: frac_origlai |
---|
1670 | |
---|
1671 | CHARACTER(LEN=80) :: meter |
---|
1672 | REAL(r_std) :: prog, sumf |
---|
1673 | LOGICAL :: found |
---|
1674 | INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog |
---|
1675 | REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max |
---|
1676 | ! |
---|
1677 | pi = 4. * ATAN(1.) |
---|
1678 | ! |
---|
1679 | !Config Key = LAI_FILE |
---|
1680 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
1681 | !Config If = !LAI_MAP |
---|
1682 | !Config Def = ../surfmap/lai2D.nc |
---|
1683 | !Config Help = The name of the file to be opened to read the LAI |
---|
1684 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
1685 | !Config map which is derived from a Nicolas VIOVY one. |
---|
1686 | ! |
---|
1687 | filename = 'lai2D.nc' |
---|
1688 | CALL getin_p('LAI_FILE',filename) |
---|
1689 | ! |
---|
1690 | IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
1691 | CALL bcast(iml) |
---|
1692 | CALL bcast(jml) |
---|
1693 | CALL bcast(lml) |
---|
1694 | CALL bcast(tml) |
---|
1695 | ! |
---|
1696 | ! |
---|
1697 | ALLOCATE(lon_lu(iml)) |
---|
1698 | ALLOCATE(lat_lu(jml)) |
---|
1699 | ALLOCATE(laimap_lu(iml,jml,nvm,tml)) |
---|
1700 | ALLOCATE(mask_lu(iml,jml)) |
---|
1701 | ! |
---|
1702 | WRITE(numout,*) 'slowproc_interlai : Reading the LAI file' |
---|
1703 | ! |
---|
1704 | IF (is_root_prc) THEN |
---|
1705 | CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu) |
---|
1706 | CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu) |
---|
1707 | CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu) |
---|
1708 | CALL flinget(fid, 'mask', iml, jml, 0, 0, 0, 1, mask_lu) |
---|
1709 | ! |
---|
1710 | CALL flinclo(fid) |
---|
1711 | ENDIF |
---|
1712 | CALL bcast(lon_lu) |
---|
1713 | CALL bcast(lat_lu) |
---|
1714 | CALL bcast(laimap_lu) |
---|
1715 | CALL bcast(mask_lu) |
---|
1716 | ! |
---|
1717 | WRITE(numout,*) 'slowproc_interlai : ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml) |
---|
1718 | ! |
---|
1719 | ! |
---|
1720 | ijml=iml*jml |
---|
1721 | ALLOCATE(lon_ful(ijml)) |
---|
1722 | ALLOCATE(lat_ful(ijml)) |
---|
1723 | ALLOCATE(laimaporig(ijml,nvm,tml)) |
---|
1724 | ALLOCATE(mask(ijml)) |
---|
1725 | |
---|
1726 | DO i=1,iml |
---|
1727 | DO j=1,jml |
---|
1728 | iki=(j-1)*iml+i |
---|
1729 | lon_ful(iki)=lon_lu(i) |
---|
1730 | lat_ful(iki)=lat_lu(j) |
---|
1731 | laimaporig(iki,:,:)=laimap_lu(i,j,:,:) |
---|
1732 | mask(iki)=mask_lu(i,j) |
---|
1733 | ENDDO |
---|
1734 | ENDDO |
---|
1735 | ! |
---|
1736 | WHERE ( laimaporig(:,:,:) .LT. 0 ) |
---|
1737 | laimaporig(:,:,:) = 0. |
---|
1738 | ENDWHERE |
---|
1739 | ! |
---|
1740 | ! |
---|
1741 | ALLOCATE(lon_up(nbpt)) |
---|
1742 | ALLOCATE(lon_low(nbpt)) |
---|
1743 | ALLOCATE(lat_up(nbpt)) |
---|
1744 | ALLOCATE(lat_low(nbpt)) |
---|
1745 | ! |
---|
1746 | DO ib =1, nbpt |
---|
1747 | ! |
---|
1748 | ! We find the 4 limits of the grid-box. As we transform the resolution of the model |
---|
1749 | ! into longitudes and latitudes we do not have the problem of periodicity. |
---|
1750 | ! coslat is a help variable here ! |
---|
1751 | ! |
---|
1752 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth |
---|
1753 | ! |
---|
1754 | lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
1755 | lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
1756 | ! |
---|
1757 | coslat = pi/180. * R_Earth |
---|
1758 | ! |
---|
1759 | lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
1760 | lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
1761 | ! |
---|
1762 | ! |
---|
1763 | ! |
---|
1764 | ENDDO |
---|
1765 | lon_up = NINT( lon_up * 1E6 ) * 1E-6 |
---|
1766 | lon_low = NINT( lon_low * 1E6 ) * 1E-6 |
---|
1767 | lat_up = NINT( lat_up * 1E6 ) * 1E-6 |
---|
1768 | lat_low = NINT( lat_low * 1E6 ) * 1E-6 |
---|
1769 | ! |
---|
1770 | ! Get the limits of the integration domaine so that we can speed up the calculations |
---|
1771 | ! |
---|
1772 | domaine_lon_min = MINVAL(lon_low) |
---|
1773 | domaine_lon_max = MAXVAL(lon_up) |
---|
1774 | domaine_lat_min = MINVAL(lat_low) |
---|
1775 | domaine_lat_max = MAXVAL(lat_up) |
---|
1776 | ! |
---|
1777 | !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max |
---|
1778 | !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max |
---|
1779 | ! |
---|
1780 | ! Ensure that the fine grid covers the whole domain |
---|
1781 | WHERE ( lon_ful(:) .LT. domaine_lon_min ) |
---|
1782 | lon_ful(:) = lon_ful(:) + 360. |
---|
1783 | ENDWHERE |
---|
1784 | ! |
---|
1785 | WHERE ( lon_ful(:) .GT. domaine_lon_max ) |
---|
1786 | lon_ful(:) = lon_ful(:) - 360. |
---|
1787 | ENDWHERE |
---|
1788 | lon_ful = NINT( lon_ful * 1E6 ) * 1E-6 |
---|
1789 | lat_ful = NINT( lat_ful * 1E6 ) * 1E-6 |
---|
1790 | ! |
---|
1791 | WRITE(numout,*) 'Interpolating the lai map :' |
---|
1792 | WRITE(numout,'(2a40)')'0%--------------------------------------', & |
---|
1793 | & '------------------------------------100%' |
---|
1794 | ! |
---|
1795 | ilast = 1 |
---|
1796 | n_origlai(:) = 0 |
---|
1797 | laimap(:,:,:) = 0. |
---|
1798 | ! |
---|
1799 | DO ip=1,ijml |
---|
1800 | ! |
---|
1801 | ! Give a progress meter |
---|
1802 | ! |
---|
1803 | ! prog = ip/float(ijml)*79. |
---|
1804 | ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(ijml)*79. ) THEN |
---|
1805 | ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' |
---|
1806 | ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) |
---|
1807 | ! WRITE(numout, advance="no", FMT='(a80)') meter |
---|
1808 | ! ENDIF |
---|
1809 | iprog = NINT(float(ip)/float(ijml)*79.) - NINT(float(ip-1)/float(ijml)*79.) |
---|
1810 | IF ( iprog .NE. 0 ) THEN |
---|
1811 | WRITE(numout,'(a1,$)') 'y' |
---|
1812 | ENDIF |
---|
1813 | ! |
---|
1814 | ! Only start looking for its place in the smaler grid if we are within the domaine |
---|
1815 | ! That should speed up things ! |
---|
1816 | ! |
---|
1817 | IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & |
---|
1818 | ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & |
---|
1819 | ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & |
---|
1820 | ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN |
---|
1821 | ! |
---|
1822 | ! look for point on GCM grid which this point on fine grid belongs to. |
---|
1823 | ! First look at the point on the model grid where we arrived just before. There is |
---|
1824 | ! a good chace that neighbouring points on the fine grid fall into the same model |
---|
1825 | ! grid box. |
---|
1826 | ! |
---|
1827 | IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. & |
---|
1828 | ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & |
---|
1829 | ( lat_ful(ip) .GE. lat_low(ilast) ) .AND. & |
---|
1830 | ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN |
---|
1831 | ! |
---|
1832 | ! We were lucky |
---|
1833 | ! |
---|
1834 | IF (mask(ip) .GT. 0) THEN |
---|
1835 | n_origlai(ilast) = n_origlai(ilast) + 1 |
---|
1836 | DO i=1,nvm |
---|
1837 | DO j=1,12 |
---|
1838 | laimap(ilast,i,j) = laimap(ilast,i,j) + laimaporig(ip,i,j) |
---|
1839 | ENDDO |
---|
1840 | ENDDO |
---|
1841 | ENDIF |
---|
1842 | ! |
---|
1843 | ELSE |
---|
1844 | ! |
---|
1845 | ! Otherwise, look everywhere. |
---|
1846 | ! Begin close to last grid point. |
---|
1847 | ! |
---|
1848 | found = .FALSE. |
---|
1849 | idi = 1 |
---|
1850 | ! |
---|
1851 | DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) |
---|
1852 | |
---|
1853 | ! |
---|
1854 | ! forward and backward |
---|
1855 | ! |
---|
1856 | DO ii = 1,2 |
---|
1857 | ! |
---|
1858 | IF ( ii .EQ. 1 ) THEN |
---|
1859 | ib = ilast - idi |
---|
1860 | ELSE |
---|
1861 | ib = ilast + idi |
---|
1862 | ENDIF |
---|
1863 | ! |
---|
1864 | IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN |
---|
1865 | IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. & |
---|
1866 | ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & |
---|
1867 | ( lat_ful(ip) .GE. lat_low(ib) ) .AND. & |
---|
1868 | ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN |
---|
1869 | ! |
---|
1870 | IF (mask(ip) .gt. 0) THEN |
---|
1871 | DO i=1,nvm |
---|
1872 | DO j=1,12 |
---|
1873 | laimap(ib,i,j) = laimap(ib,i,j) + laimaporig(ip,i,j) |
---|
1874 | ENDDO |
---|
1875 | ENDDO |
---|
1876 | n_origlai(ib) = n_origlai(ib) + 1 |
---|
1877 | ENDIF |
---|
1878 | ilast = ib |
---|
1879 | found = .TRUE. |
---|
1880 | ! |
---|
1881 | ENDIF |
---|
1882 | ENDIF |
---|
1883 | ! |
---|
1884 | ENDDO |
---|
1885 | ! |
---|
1886 | idi = idi + 1 |
---|
1887 | ! |
---|
1888 | ENDDO |
---|
1889 | ! |
---|
1890 | ENDIF ! lucky/not lucky |
---|
1891 | ! |
---|
1892 | ENDIF ! in the domain |
---|
1893 | ENDDO |
---|
1894 | |
---|
1895 | |
---|
1896 | ! determine fraction of LAI points in each box of the coarse grid |
---|
1897 | DO ip=1,nbpt |
---|
1898 | IF ( n_origlai(ip) .GT. 0 ) THEN |
---|
1899 | DO jv =1,nvm |
---|
1900 | laimap(ip,jv,:) = laimap(ip,jv,:)/REAL(n_origlai(ip),r_std) |
---|
1901 | ENDDO |
---|
1902 | ELSE |
---|
1903 | ! |
---|
1904 | ! Now we need to handle some exceptions |
---|
1905 | ! |
---|
1906 | IF ( lalo(ip,1) .LT. -56.0) THEN |
---|
1907 | ! Antartica |
---|
1908 | DO jv =1,nvm |
---|
1909 | laimap(ip,jv,:) = 0. |
---|
1910 | ENDDO |
---|
1911 | ! |
---|
1912 | ELSE IF ( lalo(ip,1) .GT. 70.0) THEN |
---|
1913 | ! Artica |
---|
1914 | DO jv =1,nvm |
---|
1915 | laimap(ip,jv,:) = 0. |
---|
1916 | ENDDO |
---|
1917 | ! |
---|
1918 | ELSE IF ( lalo(ip,1) .GT. 55.0 .AND. lalo(ip,2) .GT. -65.0 .AND. lalo(ip,2) .LT. -20.0) THEN |
---|
1919 | ! Greenland |
---|
1920 | DO jv =1,nvm |
---|
1921 | laimap(ip,jv,:) = 0. |
---|
1922 | ENDDO |
---|
1923 | ! |
---|
1924 | ELSE |
---|
1925 | ! |
---|
1926 | WRITE(numout,*) 'PROBLEM, no point in the lai map found for this grid box' |
---|
1927 | WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip) |
---|
1928 | WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip) |
---|
1929 | ! |
---|
1930 | WRITE(numout,*) 'Looking for nearest point on the lai map file' |
---|
1931 | CALL slowproc_nearest (ijml, lon_ful, lat_ful, & |
---|
1932 | lalo(ip,2), lalo(ip,1), inear) |
---|
1933 | WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', & |
---|
1934 | lon_ful(inear),lat_ful(inear) |
---|
1935 | ! |
---|
1936 | DO jv =1,nvm |
---|
1937 | laimap(ip,jv,:) = laimaporig(inear,jv,:) |
---|
1938 | ENDDO |
---|
1939 | ENDIF |
---|
1940 | ENDIF |
---|
1941 | ENDDO |
---|
1942 | ! |
---|
1943 | WRITE(numout,*) 'slowproc_interlai : Interpolation Done' |
---|
1944 | ! |
---|
1945 | ! |
---|
1946 | ! |
---|
1947 | DEALLOCATE(lon_up) |
---|
1948 | DEALLOCATE(lon_low) |
---|
1949 | DEALLOCATE(lat_up) |
---|
1950 | DEALLOCATE(lat_low) |
---|
1951 | DEALLOCATE(lat_ful) |
---|
1952 | DEALLOCATE(lon_ful) |
---|
1953 | DEALLOCATE(lat_lu) |
---|
1954 | DEALLOCATE(lon_lu) |
---|
1955 | DEALLOCATE(laimap_lu) |
---|
1956 | DEALLOCATE(laimaporig) |
---|
1957 | DEALLOCATE(mask_lu) |
---|
1958 | DEALLOCATE(mask) |
---|
1959 | ! |
---|
1960 | RETURN |
---|
1961 | ! |
---|
1962 | END SUBROUTINE slowproc_interlai_OLD |
---|
1963 | !! |
---|
1964 | !! Interpolate the LAI map to the grid of the model |
---|
1965 | !! |
---|
1966 | SUBROUTINE slowproc_interlai_NEW(nbpt, lalo, resolution, neighbours, contfrac, laimap) |
---|
1967 | ! |
---|
1968 | ! |
---|
1969 | ! |
---|
1970 | ! 0.1 INPUT |
---|
1971 | ! |
---|
1972 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
1973 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
1974 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
1975 | ! |
---|
1976 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
1977 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
1978 | ! |
---|
1979 | ! 0.2 OUTPUT |
---|
1980 | ! |
---|
1981 | REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned |
---|
1982 | ! |
---|
1983 | ! 0.3 LOCAL |
---|
1984 | ! |
---|
1985 | ! |
---|
1986 | CHARACTER(LEN=80) :: filename |
---|
1987 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv |
---|
1988 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu |
---|
1989 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon |
---|
1990 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area |
---|
1991 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index |
---|
1992 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu |
---|
1993 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask |
---|
1994 | ! |
---|
1995 | REAL(r_std) :: totarea |
---|
1996 | INTEGER(i_std) :: idi, nbvmax |
---|
1997 | CHARACTER(LEN=30) :: callsign |
---|
1998 | ! |
---|
1999 | LOGICAL :: ok_interpol ! optionnal return of aggregate_2d |
---|
2000 | ! |
---|
2001 | INTEGER :: ALLOC_ERR |
---|
2002 | ! |
---|
2003 | !Config Key = LAI_FILE |
---|
2004 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
2005 | !Config If = LAI_MAP |
---|
2006 | !Config Def = ../surfmap/lai2D.nc |
---|
2007 | !Config Help = The name of the file to be opened to read the LAI |
---|
2008 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
2009 | !Config map which is derived from a Nicolas VIOVY one. |
---|
2010 | ! |
---|
2011 | filename = 'lai2D.nc' |
---|
2012 | CALL getin_p('LAI_FILE',filename) |
---|
2013 | ! |
---|
2014 | IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
2015 | CALL bcast(iml) |
---|
2016 | CALL bcast(jml) |
---|
2017 | CALL bcast(lml) |
---|
2018 | CALL bcast(tml) |
---|
2019 | ! |
---|
2020 | ! |
---|
2021 | ALLOC_ERR=-1 |
---|
2022 | ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR) |
---|
2023 | IF (ALLOC_ERR/=0) THEN |
---|
2024 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR |
---|
2025 | STOP |
---|
2026 | ENDIF |
---|
2027 | ALLOC_ERR=-1 |
---|
2028 | ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR) |
---|
2029 | IF (ALLOC_ERR/=0) THEN |
---|
2030 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR |
---|
2031 | STOP |
---|
2032 | ENDIF |
---|
2033 | ALLOC_ERR=-1 |
---|
2034 | ALLOCATE(laimap_lu(iml,jml,nvm,tml), STAT=ALLOC_ERR) |
---|
2035 | IF (ALLOC_ERR/=0) THEN |
---|
2036 | WRITE(numout,*) "ERROR IN ALLOCATION of laimap_lu : ",ALLOC_ERR |
---|
2037 | STOP |
---|
2038 | ENDIF |
---|
2039 | ! |
---|
2040 | ! |
---|
2041 | IF (is_root_prc) THEN |
---|
2042 | CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu) |
---|
2043 | CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu) |
---|
2044 | CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu) |
---|
2045 | ! |
---|
2046 | WHERE (laimap_lu(:,:,:,:) < zero ) |
---|
2047 | laimap_lu(:,:,:,:) = zero |
---|
2048 | ENDWHERE |
---|
2049 | ! |
---|
2050 | CALL flinclo(fid) |
---|
2051 | ENDIF |
---|
2052 | CALL bcast(lon_lu) |
---|
2053 | CALL bcast(lat_lu) |
---|
2054 | CALL bcast(laimap_lu) |
---|
2055 | ! |
---|
2056 | ALLOC_ERR=-1 |
---|
2057 | ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR) |
---|
2058 | IF (ALLOC_ERR/=0) THEN |
---|
2059 | WRITE(numout,*) "ERROR IN ALLOCATION of lon : ",ALLOC_ERR |
---|
2060 | STOP |
---|
2061 | ENDIF |
---|
2062 | ALLOC_ERR=-1 |
---|
2063 | ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR) |
---|
2064 | IF (ALLOC_ERR/=0) THEN |
---|
2065 | WRITE(numout,*) "ERROR IN ALLOCATION of lat : ",ALLOC_ERR |
---|
2066 | STOP |
---|
2067 | ENDIF |
---|
2068 | ! |
---|
2069 | DO ip=1,iml |
---|
2070 | lat(ip,:) = lat_lu(:) |
---|
2071 | ENDDO |
---|
2072 | DO jp=1,jml |
---|
2073 | lon(:,jp) = lon_lu(:) |
---|
2074 | ENDDO |
---|
2075 | ! |
---|
2076 | ALLOC_ERR=-1 |
---|
2077 | ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) |
---|
2078 | IF (ALLOC_ERR/=0) THEN |
---|
2079 | WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR |
---|
2080 | STOP |
---|
2081 | ENDIF |
---|
2082 | ! |
---|
2083 | ! Consider all points a priori |
---|
2084 | ! |
---|
2085 | !!$ mask(:,:) = 1 |
---|
2086 | mask(:,:) = 0 |
---|
2087 | ! |
---|
2088 | DO ip=1,iml |
---|
2089 | DO jp=1,jml |
---|
2090 | ! |
---|
2091 | ! Exclude the points where there is never a LAI value. It is probably |
---|
2092 | ! an ocean point. |
---|
2093 | ! |
---|
2094 | !!$ IF (ABS(SUM(laimap_lu(ip,jp,:,:))) < EPSILON(laimap_lu) ) THEN |
---|
2095 | !!$ mask(ip,jp) = 0 |
---|
2096 | !!$ ENDIF |
---|
2097 | ! |
---|
2098 | IF ( ANY(laimap_lu(ip,jp,:,:) < 20.) ) THEN |
---|
2099 | mask(ip,jp) = 1 |
---|
2100 | ENDIF |
---|
2101 | ! |
---|
2102 | ENDDO |
---|
2103 | ENDDO |
---|
2104 | ! |
---|
2105 | nbvmax = 20 |
---|
2106 | ! |
---|
2107 | callsign = 'LAI map' |
---|
2108 | ! |
---|
2109 | ok_interpol = .FALSE. |
---|
2110 | DO WHILE ( .NOT. ok_interpol ) |
---|
2111 | WRITE(numout,*) "Projection arrays for ",callsign," : " |
---|
2112 | WRITE(numout,*) "nbvmax = ",nbvmax |
---|
2113 | ! |
---|
2114 | ALLOC_ERR=-1 |
---|
2115 | ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR) |
---|
2116 | IF (ALLOC_ERR/=0) THEN |
---|
2117 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR |
---|
2118 | STOP |
---|
2119 | ENDIF |
---|
2120 | sub_index(:,:,:)=0 |
---|
2121 | ALLOC_ERR=-1 |
---|
2122 | ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
2123 | IF (ALLOC_ERR/=0) THEN |
---|
2124 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR |
---|
2125 | STOP |
---|
2126 | ENDIF |
---|
2127 | sub_area(:,:)=zero |
---|
2128 | ! |
---|
2129 | CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
2130 | & iml, jml, lon, lat, mask, callsign, & |
---|
2131 | & nbvmax, sub_index, sub_area, ok_interpol) |
---|
2132 | |
---|
2133 | ! |
---|
2134 | IF ( .NOT. ok_interpol ) THEN |
---|
2135 | DEALLOCATE(sub_area) |
---|
2136 | DEALLOCATE(sub_index) |
---|
2137 | ENDIF |
---|
2138 | ! |
---|
2139 | nbvmax = nbvmax * 2 |
---|
2140 | ENDDO |
---|
2141 | ! |
---|
2142 | laimap(:,:,:) = zero |
---|
2143 | ! |
---|
2144 | DO ib=1,nbpt |
---|
2145 | idi = COUNT(sub_area(ib,:) > zero) |
---|
2146 | IF ( idi > 0 ) THEN |
---|
2147 | totarea = zero |
---|
2148 | DO jj=1,idi |
---|
2149 | ip = sub_index(ib,jj,1) |
---|
2150 | jp = sub_index(ib,jj,2) |
---|
2151 | DO jv=1,nvm |
---|
2152 | DO it=1,12 |
---|
2153 | laimap(ib,jv,it) = laimap(ib,jv,it) + laimap_lu(ip,jp,jv,it)*sub_area(ib,jj) |
---|
2154 | ENDDO |
---|
2155 | ENDDO |
---|
2156 | totarea = totarea + sub_area(ib,jj) |
---|
2157 | ENDDO |
---|
2158 | ! |
---|
2159 | ! Normalize |
---|
2160 | ! |
---|
2161 | laimap(ib,:,:) = laimap(ib,:,:)/totarea |
---|
2162 | |
---|
2163 | !!$ IF ( ANY( laimap(ib,:,:) > 10 ) ) THEN |
---|
2164 | !!$ WRITE(numout,*) "For point ",ib |
---|
2165 | !!$ WRITE(numout,*) lalo(ib,1), lalo(ib,2) |
---|
2166 | !!$ WRITE(numout,*) "For ib=",ib," we have LAI ",laimap(ib,:,1:idi) |
---|
2167 | !!$ WRITE(numout,*) "Detail of projection :" |
---|
2168 | !!$ WRITE(numout,*) sub_area(ib,1:idi) |
---|
2169 | !!$ WRITE(numout,*) "Total for projection :" ,totarea |
---|
2170 | !!$ ENDIF |
---|
2171 | ! |
---|
2172 | ELSE |
---|
2173 | WRITE(numout,*) 'On point ', ib, ' no points where found for interpolating the LAI map.' |
---|
2174 | WRITE(numout,*) 'Location : ', lalo(ib,2), lalo(ib,1) |
---|
2175 | DO jv=1,nvm |
---|
2176 | laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux |
---|
2177 | ENDDO |
---|
2178 | WRITE(numout,*) 'Solved by putting the average LAI for the PFT all year long' |
---|
2179 | ENDIF |
---|
2180 | ENDDO |
---|
2181 | ! |
---|
2182 | WRITE(numout,*) 'slowproc_interlai : Interpolation Done' |
---|
2183 | ! |
---|
2184 | ! |
---|
2185 | ! |
---|
2186 | DEALLOCATE(lat_lu) |
---|
2187 | DEALLOCATE(lon_lu) |
---|
2188 | DEALLOCATE(lon) |
---|
2189 | DEALLOCATE(lat) |
---|
2190 | DEALLOCATE(laimap_lu) |
---|
2191 | DEALLOCATE(mask) |
---|
2192 | DEALLOCATE(sub_area) |
---|
2193 | DEALLOCATE(sub_index) |
---|
2194 | ! |
---|
2195 | RETURN |
---|
2196 | ! |
---|
2197 | END SUBROUTINE slowproc_interlai_NEW |
---|
2198 | !! |
---|
2199 | !! Interpolate a vegetation map (by pft) |
---|
2200 | !! |
---|
2201 | !MM modif TAG 1.4 : |
---|
2202 | ! SUBROUTINE slowproc_update(nbpt, lalo, resolution, vegetmap, frac_nobiomap) |
---|
2203 | !MM modif TAG 1.9.2 : |
---|
2204 | ! SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, vegetmax, frac_nobio, veget_year, init) |
---|
2205 | SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
2206 | & veget_last, frac_nobio_last, & |
---|
2207 | & veget_next, frac_nobio_next, veget_year, init) |
---|
2208 | ! |
---|
2209 | ! |
---|
2210 | ! |
---|
2211 | ! 0.1 INPUT |
---|
2212 | ! |
---|
2213 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs |
---|
2214 | ! to be interpolated |
---|
2215 | REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo ! Vector of latitude and longitudes (beware of the order !) |
---|
2216 | !MM modif TAG 1.4 : add grid variables for aggregate |
---|
2217 | INTEGER(i_std), DIMENSION(nbpt,8), INTENT(in) :: neighbours ! Vector of neighbours for each grid point |
---|
2218 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
2219 | REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution ! The size in km of each grid-box in X and Y |
---|
2220 | REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of continent in the grid |
---|
2221 | ! |
---|
2222 | REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last ! old max vegetfrac |
---|
2223 | REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(in) :: frac_nobio_last ! old fraction of the mesh which is |
---|
2224 | ! covered by ice, lakes, ... |
---|
2225 | ! |
---|
2226 | INTEGER(i_std), INTENT(in) :: veget_year ! first year for landuse (0 == NO TIME AXIS) |
---|
2227 | LOGICAL, OPTIONAL, INTENT(in) :: init ! initialisation : in case of dgvm, it forces update of all PFTs |
---|
2228 | ! |
---|
2229 | ! 0.2 OUTPUT |
---|
2230 | ! |
---|
2231 | REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next ! new max vegetfrac |
---|
2232 | REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next ! new fraction of the mesh which is |
---|
2233 | ! covered by ice, lakes, ... |
---|
2234 | ! |
---|
2235 | ! 0.3 LOCAL |
---|
2236 | ! |
---|
2237 | ! |
---|
2238 | CHARACTER(LEN=80) :: filename |
---|
2239 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, inobio, jv |
---|
2240 | INTEGER(i_std) :: nb_coord,nb_var, nb_gat,nb_dim |
---|
2241 | REAL(r_std) :: date, dt |
---|
2242 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: itau |
---|
2243 | REAL(r_std), DIMENSION(1) :: time_counter |
---|
2244 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu |
---|
2245 | INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w, i_d_w |
---|
2246 | LOGICAL :: exv, l_ex |
---|
2247 | ! |
---|
2248 | !MM modif TAG 1.4 : suppression of all time axis reading and interpolation, replaced by each year 2D reading. |
---|
2249 | ! REAL(r_std), INTENT(inout) :: vegetmap(nbpt,nvm,293) ! vegetfrac read variable and re-dimensioned |
---|
2250 | ! REAL(r_std), INTENT(inout) :: frac_nobiomap(nbpt,nnobio,293) ! Fraction of the mesh which is covered by ice, lakes, ... |
---|
2251 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: vegmap ! last coord is time with only one value = 1 |
---|
2252 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: vegmap_1 ! last coord is time with only one value = 1 (IF VEGET_YEAR == 0 , NO TIME AXIS) |
---|
2253 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful |
---|
2254 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area |
---|
2255 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index |
---|
2256 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask |
---|
2257 | ! |
---|
2258 | REAL(r_std) :: sumf, err, norm |
---|
2259 | REAL(r_std) :: totarea |
---|
2260 | INTEGER(i_std) :: idi, nbvmax |
---|
2261 | CHARACTER(LEN=30) :: callsign |
---|
2262 | ! |
---|
2263 | LOGICAL :: ok_interpol ! optionnal return of aggregate_2d |
---|
2264 | ! |
---|
2265 | ! for DGVM case : |
---|
2266 | REAL(r_std) :: sum_veg ! sum of vegets |
---|
2267 | REAL(r_std) :: sum_nobio ! sum of nobios |
---|
2268 | REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets |
---|
2269 | REAL(r_std) :: rapport ! (S-B) / (S-A) |
---|
2270 | LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones) |
---|
2271 | ! e.g. in case of DGVM and not init (optional parameter) |
---|
2272 | ! |
---|
2273 | LOGICAL, PARAMETER :: debug = .FALSE. |
---|
2274 | ! |
---|
2275 | INTEGER :: ALLOC_ERR |
---|
2276 | ! |
---|
2277 | !Config Key = VEGETATION_FILE |
---|
2278 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
2279 | !Config If = LAND_USE |
---|
2280 | !Config Def = pft_new.nc |
---|
2281 | !Config Help = The name of the file to be opened to read a vegetation |
---|
2282 | !Config map (in pft) is to be given here. |
---|
2283 | ! |
---|
2284 | filename = 'pft_new.nc' |
---|
2285 | CALL getin_p('VEGETATION_FILE',filename) |
---|
2286 | ! |
---|
2287 | IF (is_root_prc) THEN |
---|
2288 | IF (debug) THEN |
---|
2289 | WRITE(numout,*) "Entering slowproc_update. Debug mode." |
---|
2290 | WRITE (*,'(/," --> fliodmpf")') |
---|
2291 | CALL fliodmpf (TRIM(filename)) |
---|
2292 | WRITE (*,'(/," --> flioopfd")') |
---|
2293 | ENDIF |
---|
2294 | CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat) |
---|
2295 | IF (debug) THEN |
---|
2296 | WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord |
---|
2297 | WRITE (*,'(" Number of variables in the file : ",I2)') nb_var |
---|
2298 | WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat |
---|
2299 | ENDIF |
---|
2300 | ENDIF |
---|
2301 | CALL bcast(nb_coord) |
---|
2302 | CALL bcast(nb_var) |
---|
2303 | CALL bcast(nb_gat) |
---|
2304 | IF ( veget_year > 0 ) THEN |
---|
2305 | IF (is_root_prc) & |
---|
2306 | CALL flioinqv (fid,v_n="time_counter",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) |
---|
2307 | CALL bcast(l_d_w) |
---|
2308 | tml=l_d_w(1) |
---|
2309 | WRITE(numout,*) "tml =",tml |
---|
2310 | ENDIF |
---|
2311 | IF (is_root_prc) & |
---|
2312 | CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) |
---|
2313 | CALL bcast(l_d_w) |
---|
2314 | iml=l_d_w(1) |
---|
2315 | WRITE(numout,*) "iml =",iml |
---|
2316 | |
---|
2317 | IF (is_root_prc) & |
---|
2318 | CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) |
---|
2319 | CALL bcast(l_d_w) |
---|
2320 | jml=l_d_w(1) |
---|
2321 | WRITE(numout,*) "jml =",jml |
---|
2322 | |
---|
2323 | IF (is_root_prc) & |
---|
2324 | CALL flioinqv (fid,v_n="veget",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) |
---|
2325 | CALL bcast(l_d_w) |
---|
2326 | lml=l_d_w(1) |
---|
2327 | |
---|
2328 | IF (lml /= nvm) & |
---|
2329 | CALL ipslerr (3,'slowproc_update', & |
---|
2330 | & 'Problem with vegetation file for Land Use','lml /= nvm', & |
---|
2331 | & '(number of pft must be equal)') |
---|
2332 | ! |
---|
2333 | ALLOC_ERR=-1 |
---|
2334 | ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR) |
---|
2335 | IF (ALLOC_ERR/=0) THEN |
---|
2336 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR |
---|
2337 | STOP |
---|
2338 | ENDIF |
---|
2339 | ALLOC_ERR=-1 |
---|
2340 | ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR) |
---|
2341 | IF (ALLOC_ERR/=0) THEN |
---|
2342 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR |
---|
2343 | STOP |
---|
2344 | ENDIF |
---|
2345 | |
---|
2346 | IF ( veget_year > 0 ) THEN |
---|
2347 | IF (tml > 0) THEN |
---|
2348 | ALLOC_ERR=-1 |
---|
2349 | ALLOCATE(itau(tml), STAT=ALLOC_ERR) |
---|
2350 | IF (ALLOC_ERR/=0) THEN |
---|
2351 | WRITE(numout,*) "ERROR IN ALLOCATION of itau : ",ALLOC_ERR |
---|
2352 | STOP |
---|
2353 | ENDIF |
---|
2354 | ENDIF |
---|
2355 | ! |
---|
2356 | IF (is_root_prc) THEN |
---|
2357 | IF (tml > 0) THEN |
---|
2358 | CALL fliogstc (fid, t_axis=itau,x_axis=lon_lu,y_axis=lat_lu) |
---|
2359 | IF (veget_year <= tml) THEN |
---|
2360 | CALL fliogetv (fid,"time_counter",time_counter, start=(/ veget_year /), count=(/ 1 /)) |
---|
2361 | WRITE(numout,*) "slowproc_update LAND_USE : the date read for vegetmax is ",time_counter |
---|
2362 | ELSE |
---|
2363 | CALL fliogetv (fid,"time_counter",time_counter, start=(/ tml /), count=(/ 1 /)) |
---|
2364 | WRITE(numout,*) "slowproc_update LAND_USE : You try to update vegetmax with a the date greater than in the file." |
---|
2365 | WRITE(numout,*) " We will keep the last one :",time_counter |
---|
2366 | ENDIF |
---|
2367 | ELSE |
---|
2368 | CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu) |
---|
2369 | ENDIF |
---|
2370 | ENDIF |
---|
2371 | ENDIF |
---|
2372 | IF (tml > 0) THEN |
---|
2373 | CALL bcast(itau) |
---|
2374 | ENDIF |
---|
2375 | CALL bcast(lon_lu) |
---|
2376 | CALL bcast(lat_lu) |
---|
2377 | ! |
---|
2378 | ALLOC_ERR=-1 |
---|
2379 | ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR) |
---|
2380 | IF (ALLOC_ERR/=0) THEN |
---|
2381 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR |
---|
2382 | STOP |
---|
2383 | ENDIF |
---|
2384 | ALLOC_ERR=-1 |
---|
2385 | ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR) |
---|
2386 | IF (ALLOC_ERR/=0) THEN |
---|
2387 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR |
---|
2388 | STOP |
---|
2389 | ENDIF |
---|
2390 | ! |
---|
2391 | DO ip=1,iml |
---|
2392 | lon_ful(ip,:)=lon_lu(ip) |
---|
2393 | ENDDO |
---|
2394 | DO jp=1,jml |
---|
2395 | lat_ful(:,jp)=lat_lu(jp) |
---|
2396 | ENDDO |
---|
2397 | ! |
---|
2398 | ! |
---|
2399 | WRITE(numout,*) 'Reading the LAND USE vegetation file' |
---|
2400 | ! |
---|
2401 | ALLOC_ERR=-1 |
---|
2402 | ALLOCATE(vegmap(iml,jml,nvm,1), STAT=ALLOC_ERR) |
---|
2403 | IF (ALLOC_ERR/=0) THEN |
---|
2404 | WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR |
---|
2405 | STOP |
---|
2406 | ENDIF |
---|
2407 | IF ( veget_year == 0 ) THEN |
---|
2408 | IF (is_root_prc) THEN |
---|
2409 | ALLOC_ERR=-1 |
---|
2410 | ALLOCATE(vegmap_1(iml,jml,nvm), STAT=ALLOC_ERR) |
---|
2411 | IF (ALLOC_ERR/=0) THEN |
---|
2412 | WRITE(numout,*) "ERROR IN ALLOCATION of vegmap_1 : ",ALLOC_ERR |
---|
2413 | STOP |
---|
2414 | ENDIF |
---|
2415 | ENDIF |
---|
2416 | ENDIF |
---|
2417 | ! |
---|
2418 | !!$ CALL flinopen & |
---|
2419 | !!$ & (filename, .FALSE., iml, jml, lml, lon_ful, lat_ful, & |
---|
2420 | !!$ & lev_ful, tml, itau, date, dt, fid) |
---|
2421 | !=> FATAL ERROR FROM ROUTINE flinopen |
---|
2422 | ! --> No time axis found |
---|
2423 | |
---|
2424 | !MM modif TAG 1.4 : |
---|
2425 | ! CALL flinget(fid, 'lon', iml, 0, 0, 0, 1, 1, lon_lu) |
---|
2426 | ! CALL flinget(fid, 'lat', jml, 0, 0, 0, 1, 1, lat_lu) |
---|
2427 | ! CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, 1, 293, vegmap_lu) |
---|
2428 | !FATAL ERROR FROM ROUTINE flinopen |
---|
2429 | ! --> No variable lon |
---|
2430 | ! We get only the right year |
---|
2431 | !!$ CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, veget_year, veget_year, vegmap) |
---|
2432 | !!$ ! |
---|
2433 | !!$ CALL flinclo(fid) |
---|
2434 | |
---|
2435 | IF (is_root_prc) & |
---|
2436 | CALL flioinqv (fid,"maxvegetfrac",exv,nb_dims=nb_dim,len_dims=l_d_w,id_dims=i_d_w) |
---|
2437 | CALL bcast(exv) |
---|
2438 | CALL bcast(nb_dim) |
---|
2439 | CALL bcast(l_d_w) |
---|
2440 | CALL bcast(i_d_w) |
---|
2441 | |
---|
2442 | IF (exv) THEN |
---|
2443 | IF (debug) THEN |
---|
2444 | WRITE (numout,'(" Number of dimensions : ",I2)') nb_dim |
---|
2445 | WRITE (numout,'(" Dimensions :",/,5(1X,I7,:))') l_d_w(1:nb_dim) |
---|
2446 | WRITE (numout,'(" Identifiers :",/,5(1X,I7,:))') i_d_w(1:nb_dim) |
---|
2447 | ENDIF |
---|
2448 | ! |
---|
2449 | IF ( veget_year > 0 ) THEN |
---|
2450 | IF (is_root_prc) THEN |
---|
2451 | IF (veget_year <= tml) THEN |
---|
2452 | CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, veget_year /), count=(/ iml, jml, nvm, 1 /)) |
---|
2453 | ELSE |
---|
2454 | CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, tml /), count=(/ iml, jml, nvm, 1 /)) |
---|
2455 | ENDIF |
---|
2456 | ENDIF |
---|
2457 | ELSE |
---|
2458 | IF (is_root_prc) THEN |
---|
2459 | CALL fliogetv (fid,"maxvegetfrac", vegmap_1, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvm /)) |
---|
2460 | vegmap(:,:,:,1)=vegmap_1(:,:,:) |
---|
2461 | DEALLOCATE(vegmap_1) |
---|
2462 | ENDIF |
---|
2463 | ENDIF |
---|
2464 | CALL bcast(vegmap) |
---|
2465 | IF (is_root_prc) CALL flioclo (fid) |
---|
2466 | ELSE |
---|
2467 | CALL ipslerr (3,'slowproc_update', & |
---|
2468 | & 'Problem with vegetation file for Land Use.', & |
---|
2469 | & "Variable maxvegetfrac doesn't exist.", & |
---|
2470 | & '(verify your land use file.)') |
---|
2471 | ENDIF |
---|
2472 | ! |
---|
2473 | ! Mask of permitted variables. |
---|
2474 | ! |
---|
2475 | ALLOC_ERR=-1 |
---|
2476 | ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) |
---|
2477 | IF (ALLOC_ERR/=0) THEN |
---|
2478 | WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR |
---|
2479 | STOP |
---|
2480 | ENDIF |
---|
2481 | ! |
---|
2482 | mask(:,:) = 0 |
---|
2483 | DO ip=1,iml |
---|
2484 | DO jp=1,jml |
---|
2485 | sum_veg=SUM(vegmap(ip,jp,:,1)) |
---|
2486 | IF ( sum_veg .GE. min_sechiba .AND. sum_veg .LE. 1.-1.e-7) THEN |
---|
2487 | mask(ip,jp) = 1 |
---|
2488 | IF (debug) THEN |
---|
2489 | WRITE(numout,*) "update : SUM(vegmap(",ip,jp,")) = ",sum_veg |
---|
2490 | ENDIF |
---|
2491 | ELSEIF ( sum_veg .GT. 1.-1.e-7 .AND. sum_veg .LE. 2.) THEN |
---|
2492 | ! normalization |
---|
2493 | vegmap(ip,jp,:,1) = vegmap(ip,jp,:,1) / sum_veg |
---|
2494 | mask(ip,jp) = 1 |
---|
2495 | IF (debug) THEN |
---|
2496 | WRITE(numout,*) "update : SUM(vegmap(",ip,jp,"))_c = ",SUM(vegmap(ip,jp,:,1)) |
---|
2497 | ENDIF |
---|
2498 | ENDIF |
---|
2499 | ENDDO |
---|
2500 | ENDDO |
---|
2501 | ! |
---|
2502 | ! |
---|
2503 | ! The number of maximum vegetation map points in the GCM grid should |
---|
2504 | ! also be computed and not imposed here. |
---|
2505 | ! |
---|
2506 | nbvmax = 200 |
---|
2507 | ! |
---|
2508 | callsign="Land Use Vegetation map" |
---|
2509 | ! |
---|
2510 | ok_interpol = .FALSE. |
---|
2511 | DO WHILE ( .NOT. ok_interpol ) |
---|
2512 | WRITE(numout,*) "Projection arrays for ",callsign," : " |
---|
2513 | WRITE(numout,*) "nbvmax = ",nbvmax |
---|
2514 | ! |
---|
2515 | ALLOC_ERR=-1 |
---|
2516 | ALLOCATE(sub_index(nbpt, nbvmax,2), STAT=ALLOC_ERR) |
---|
2517 | IF (ALLOC_ERR/=0) THEN |
---|
2518 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR |
---|
2519 | STOP |
---|
2520 | ENDIF |
---|
2521 | sub_index(:,:,:)=0 |
---|
2522 | |
---|
2523 | ALLOC_ERR=-1 |
---|
2524 | ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
2525 | IF (ALLOC_ERR/=0) THEN |
---|
2526 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR |
---|
2527 | STOP |
---|
2528 | ENDIF |
---|
2529 | sub_area(:,:)=zero |
---|
2530 | ! |
---|
2531 | CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
2532 | & iml, jml, lon_ful, lat_ful, mask, callsign, & |
---|
2533 | & nbvmax, sub_index, sub_area, ok_interpol) |
---|
2534 | ! |
---|
2535 | IF ( .NOT. ok_interpol ) THEN |
---|
2536 | DEALLOCATE(sub_area) |
---|
2537 | DEALLOCATE(sub_index) |
---|
2538 | ENDIF |
---|
2539 | ! |
---|
2540 | nbvmax = nbvmax * 2 |
---|
2541 | ENDDO |
---|
2542 | ! |
---|
2543 | ! Compute the logical for partial (only anthropic) PTFs update |
---|
2544 | IF (PRESENT(init)) THEN |
---|
2545 | partial_update = control%ok_dgvm .AND. .NOT. init |
---|
2546 | ELSE |
---|
2547 | partial_update = control%ok_dgvm |
---|
2548 | ENDIF |
---|
2549 | ! |
---|
2550 | IF ( .NOT. partial_update ) THEN |
---|
2551 | ! |
---|
2552 | veget_next(:,:)=zero |
---|
2553 | ! |
---|
2554 | DO ib = 1, nbpt |
---|
2555 | idi=1 |
---|
2556 | sumf=0. |
---|
2557 | DO WHILE ( sub_area(ib,idi) > zero ) |
---|
2558 | ip = sub_index(ib,idi,1) |
---|
2559 | jp = sub_index(ib,idi,2) |
---|
2560 | veget_next(ib,:) = veget_next(ib,:) + vegmap(ip,jp,:,1)*sub_area(ib,idi) |
---|
2561 | sumf=sumf + sub_area(ib,idi) |
---|
2562 | idi = idi +1 |
---|
2563 | ENDDO |
---|
2564 | !!$ ! |
---|
2565 | !!$ ! Limit the smalest vegetation fraction to 0.5% |
---|
2566 | !!$ ! |
---|
2567 | !!$ DO jv = 1, nvm |
---|
2568 | !!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN |
---|
2569 | !!$ veget_next(ib,jv) = zero |
---|
2570 | !!$ ENDIF |
---|
2571 | !!$ ENDDO |
---|
2572 | ! |
---|
2573 | ! Normalize |
---|
2574 | ! |
---|
2575 | IF (sumf > min_sechiba) THEN |
---|
2576 | veget_next(ib,:) = veget_next(ib,:) / sumf |
---|
2577 | ELSE |
---|
2578 | WRITE(numout,*) "slowproc_update : No land point in the map for point ",& |
---|
2579 | ib,",(",lalo(ib,1),",",lalo(ib,2),")" |
---|
2580 | CALL ipslerr (2,'slowproc_update', & |
---|
2581 | & 'Problem with vegetation file for Land Use.', & |
---|
2582 | & "No land point in the map for point", & |
---|
2583 | & 'Keep old values. (verify your land use file.)') |
---|
2584 | !!$ CALL slowproc_nearest (iml, lon_ful, lat_ful, & |
---|
2585 | !!$ lalo(ib,2), lalo(ib,1), inear) |
---|
2586 | IF (PRESENT(init)) THEN |
---|
2587 | IF (init) THEN |
---|
2588 | veget_next(ib,:) = (/ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. /) |
---|
2589 | ELSE |
---|
2590 | veget_next(ib,:) = veget_last(ib,:) |
---|
2591 | ENDIF |
---|
2592 | ELSE |
---|
2593 | veget_next(ib,:) = veget_last(ib,:) |
---|
2594 | ENDIF |
---|
2595 | ENDIF |
---|
2596 | ! |
---|
2597 | IF (debug) THEN |
---|
2598 | WRITE(numout,*) "SUM(veget_next(",ib,")) = ",SUM(veget_next(ib,:)) |
---|
2599 | ENDIF |
---|
2600 | ENDDO |
---|
2601 | ELSE |
---|
2602 | DO ib = 1, nbpt |
---|
2603 | ! last veget for this point |
---|
2604 | sum_veg=SUM(veget_last(ib,:)) |
---|
2605 | IF (debug) THEN |
---|
2606 | WRITE(*,*) "SUM(veget_last(",ib,")) = ",sum_veg |
---|
2607 | ENDIF |
---|
2608 | ! |
---|
2609 | ! If the DGVM is activated, only anthropiques PFT are utpdated, |
---|
2610 | ! other are copied |
---|
2611 | veget_next(ib,:) = veget_last(ib,:) |
---|
2612 | ! |
---|
2613 | ! natural ones are initialized to zero. |
---|
2614 | DO jv = 2, nvm |
---|
2615 | ! If the DGVM is activated, only anthropiques PFT are utpdated |
---|
2616 | IF ( .NOT. natural(jv) ) THEN |
---|
2617 | veget_next(ib,jv) = zero |
---|
2618 | ENDIF |
---|
2619 | ENDDO |
---|
2620 | ! |
---|
2621 | idi=1 |
---|
2622 | sumf=0. |
---|
2623 | DO WHILE ( sub_area(ib,idi) > zero ) |
---|
2624 | ip = sub_index(ib,idi,1) |
---|
2625 | jp = sub_index(ib,idi,2) |
---|
2626 | ! If the DGVM is activated, only anthropic PFTs are utpdated |
---|
2627 | DO jv = 2, nvm |
---|
2628 | IF ( .NOT. natural(jv) ) THEN |
---|
2629 | veget_next(ib,jv) = veget_next(ib,jv) + vegmap(ip,jp,jv,1)*sub_area(ib,idi) |
---|
2630 | ENDIF |
---|
2631 | ENDDO |
---|
2632 | sumf=sumf + sub_area(ib,idi) |
---|
2633 | idi = idi +1 |
---|
2634 | ENDDO |
---|
2635 | !!$ ! |
---|
2636 | !!$ ! Limit the smalest vegetation fraction to 0.5% |
---|
2637 | !!$ ! |
---|
2638 | !!$ DO jv = 2, nvm |
---|
2639 | !!$ ! On anthropic and natural PFTs ? |
---|
2640 | !!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN |
---|
2641 | !!$ veget_next(ib,jv) = zero |
---|
2642 | !!$ ENDIF |
---|
2643 | !!$ ENDDO |
---|
2644 | ! |
---|
2645 | ! Normalize |
---|
2646 | ! |
---|
2647 | ! Proposition de Pierre : |
---|
2648 | ! apres modification de la surface des PFTs anthropiques, |
---|
2649 | ! on doit conserver la proportion des PFTs naturels. |
---|
2650 | ! ie la somme des vegets est conservee |
---|
2651 | ! et PFT naturel / (somme des vegets - somme des vegets anthropiques) |
---|
2652 | ! est conservee. |
---|
2653 | ! Sum veget_next = old (sum veget_next Naturel) + (sum veget_next Anthropic) |
---|
2654 | ! = new (sum veget_next Naturel) + (sum veget_next Anthropic) |
---|
2655 | ! a / (S-A) = e / (S-B) ; b/(S-A) = f/(S-B) |
---|
2656 | IF (sumf > min_sechiba) THEN |
---|
2657 | sumvAnthro_old = zero |
---|
2658 | sumvAnthro = zero |
---|
2659 | DO jv = 2, nvm |
---|
2660 | IF ( .NOT. natural(jv) ) THEN |
---|
2661 | veget_next(ib,jv) = veget_next(ib,jv) / sumf |
---|
2662 | sumvAnthro = sumvAnthro + veget_last(ib,jv) |
---|
2663 | sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv) |
---|
2664 | ENDIF |
---|
2665 | ENDDO |
---|
2666 | ! conservation : |
---|
2667 | rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old ) |
---|
2668 | veget_next(ib,1) = veget_last(ib,1) * rapport |
---|
2669 | DO jv = 2, nvm |
---|
2670 | IF ( .NOT. natural(jv) ) THEN |
---|
2671 | veget_next(ib,jv) = veget_last(ib,jv) * rapport |
---|
2672 | ENDIF |
---|
2673 | ENDDO |
---|
2674 | ! test |
---|
2675 | IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > EPSILON(un) ) THEN |
---|
2676 | WRITE(numout,*) "No conservation of sum of veget for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" |
---|
2677 | WRITE(numout,*) "last sum of veget ",sum_veg," new sum of veget ",SUM(veget_next(ib,:))," error : ",& |
---|
2678 | & SUM(veget_next(ib,:)) - sum_veg |
---|
2679 | WRITE(numout,*) "Anthropic modificaztions : last ",sumvAnthro_old," new ",sumvAnthro |
---|
2680 | CALL ipslerr (3,'slowproc_update', & |
---|
2681 | & 'No conservation of sum of veget_next', & |
---|
2682 | & "The sum of veget_next is different after reading Land Use map.", & |
---|
2683 | & '(verify the dgvm case model.)') |
---|
2684 | ENDIF |
---|
2685 | ELSE |
---|
2686 | WRITE(numout,*) "No land point in the map for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" |
---|
2687 | ! CALL ipslerr (3,'slowproc_update', & |
---|
2688 | CALL ipslerr (2,'slowproc_update', & |
---|
2689 | & 'Problem with vegetation file for Land Use.', & |
---|
2690 | & "No land point in the map for point", & |
---|
2691 | & '(verify your land use file.)') |
---|
2692 | veget_next(ib,:) = veget_last(ib,:) |
---|
2693 | ENDIF |
---|
2694 | |
---|
2695 | ENDDO |
---|
2696 | ENDIF |
---|
2697 | ! |
---|
2698 | frac_nobio_next (:,:) = un |
---|
2699 | ! |
---|
2700 | !MM |
---|
2701 | ! Work only for one nnobio !! (ie ice) |
---|
2702 | DO inobio=1,nnobio |
---|
2703 | DO jv=1,nvm |
---|
2704 | ! |
---|
2705 | DO ib = 1, nbpt |
---|
2706 | frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv) |
---|
2707 | ENDDO |
---|
2708 | ENDDO |
---|
2709 | ! |
---|
2710 | ENDDO |
---|
2711 | ! |
---|
2712 | DO ib = 1, nbpt |
---|
2713 | sum_veg = SUM(veget_next(ib,:)) |
---|
2714 | sum_nobio = SUM(frac_nobio_next(ib,:)) |
---|
2715 | IF (sum_nobio < 0.) THEN |
---|
2716 | frac_nobio_next(ib,:) = zero |
---|
2717 | veget_next(ib,1) = veget_next(ib,1) - sum_nobio |
---|
2718 | sum_veg = SUM(veget_next(ib,:)) |
---|
2719 | ENDIF |
---|
2720 | sumf = sum_veg + sum_nobio |
---|
2721 | IF (sumf > min_sechiba) THEN |
---|
2722 | veget_next(ib,:) = veget_next(ib,:) / sumf |
---|
2723 | frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf |
---|
2724 | norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:)) |
---|
2725 | err=norm-un |
---|
2726 | IF (debug) & |
---|
2727 | WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-1., sumf",err,sumf |
---|
2728 | IF (abs(err) > -EPSILON(1._r_std)) THEN |
---|
2729 | !MM 1.9.3 |
---|
2730 | ! IF (abs(err) > 0.) THEN |
---|
2731 | IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN |
---|
2732 | frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err |
---|
2733 | ELSE |
---|
2734 | veget_next(ib,1) = veget_next(ib,1) - err |
---|
2735 | ENDIF |
---|
2736 | norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:)) |
---|
2737 | err=norm-un |
---|
2738 | IF (debug) & |
---|
2739 | WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-1.",err |
---|
2740 | IF (abs(err) > EPSILON(1._r_std)) THEN |
---|
2741 | !MM 1.9.3 |
---|
2742 | ! IF (abs(err) > 0.) THEN |
---|
2743 | WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" |
---|
2744 | WRITE(numout,*) " err(sum-1.) = ",abs(err) |
---|
2745 | CALL ipslerr (2,'slowproc_update', & |
---|
2746 | & 'Problem with sum vegetation + sum fracnobio for Land Use.', & |
---|
2747 | & "sum not equal to 1.", & |
---|
2748 | & '(verify your land use file.)') |
---|
2749 | ENDIF |
---|
2750 | ENDIF |
---|
2751 | ELSE |
---|
2752 | WRITE(numout,*) "No vegetation nor frac_nobio for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" |
---|
2753 | WRITE(numout,*)"Replaced by bare_soil !! " |
---|
2754 | veget_next(ib,1) = un |
---|
2755 | veget_next(ib,2:nvm) = zero |
---|
2756 | frac_nobio_next(ib,:) = zero |
---|
2757 | !!$ CALL ipslerr (3,'slowproc_update', & |
---|
2758 | !!$ & 'Problem with vegetation file for Land Use.', & |
---|
2759 | !!$ & "No vegetation nor frac_nobio for point ", & |
---|
2760 | !!$ & '(verify your land use file.)') |
---|
2761 | ENDIF |
---|
2762 | ENDDO |
---|
2763 | ! |
---|
2764 | WRITE(numout,*) 'slowproc_update : Interpolation Done' |
---|
2765 | ! |
---|
2766 | DEALLOCATE(vegmap) |
---|
2767 | DEALLOCATE(lat_lu,lon_lu) |
---|
2768 | DEALLOCATE(lat_ful,lon_ful) |
---|
2769 | DEALLOCATE(mask) |
---|
2770 | DEALLOCATE(sub_index,sub_area) |
---|
2771 | ! |
---|
2772 | RETURN |
---|
2773 | ! |
---|
2774 | END SUBROUTINE slowproc_update |
---|
2775 | |
---|
2776 | !! |
---|
2777 | !! Interpolate the IGBP vegetation map to the grid of the model |
---|
2778 | !MM TAG 1.6 model ! |
---|
2779 | !! |
---|
2780 | SUBROUTINE slowproc_interpol_OLD(nbpt, lalo, neighbours, resolution, veget, frac_nobio ) |
---|
2781 | ! |
---|
2782 | ! |
---|
2783 | ! |
---|
2784 | ! 0.1 INPUT |
---|
2785 | ! |
---|
2786 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
2787 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
2788 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point |
---|
2789 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
2790 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
2791 | ! |
---|
2792 | ! 0.2 OUTPUT |
---|
2793 | ! |
---|
2794 | REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions |
---|
2795 | REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... |
---|
2796 | ! |
---|
2797 | ! 0.3 LOCAL |
---|
2798 | ! |
---|
2799 | REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8 |
---|
2800 | INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes |
---|
2801 | ! |
---|
2802 | ! |
---|
2803 | CHARACTER(LEN=80) :: filename |
---|
2804 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid |
---|
2805 | REAL(r_std) :: lev(1), date, dt, coslat, pi |
---|
2806 | INTEGER(i_std) :: itau(1) |
---|
2807 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap |
---|
2808 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low |
---|
2809 | INTEGER, DIMENSION(nbpt,nolson) :: n_origveg |
---|
2810 | INTEGER, DIMENSION(nbpt) :: n_found |
---|
2811 | REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg |
---|
2812 | REAL(r_std) :: vegcorr(nolson,nvm) |
---|
2813 | REAL(r_std) :: nobiocorr(nolson,nnobio) |
---|
2814 | CHARACTER(LEN=80) :: meter |
---|
2815 | REAL(r_std) :: prog, sumf |
---|
2816 | LOGICAL :: found |
---|
2817 | INTEGER :: idi, ilast, ii, jv, inear, iprog |
---|
2818 | REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max |
---|
2819 | ! |
---|
2820 | pi = 4. * ATAN(1.) |
---|
2821 | ! |
---|
2822 | CALL get_vegcorr (nolson,vegcorr,nobiocorr) |
---|
2823 | ! |
---|
2824 | !Config Key = VEGETATION_FILE |
---|
2825 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
2826 | !Config If = !IMPOSE_VEG |
---|
2827 | !Config Def = ../surfmap/carteveg5km.nc |
---|
2828 | !Config Help = The name of the file to be opened to read the vegetation |
---|
2829 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
2830 | !Config map which is derived from the IGBP one. We assume that we have |
---|
2831 | !Config a classification in 87 types. This is Olson modified by Viovy. |
---|
2832 | ! |
---|
2833 | filename = '../surfmap/carteveg5km.nc' |
---|
2834 | CALL getin_p('VEGETATION_FILE',filename) |
---|
2835 | ! |
---|
2836 | if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
2837 | CALL bcast(iml) |
---|
2838 | CALL bcast(jml) |
---|
2839 | CALL bcast(lml) |
---|
2840 | CALL bcast(tml) |
---|
2841 | ! |
---|
2842 | ! |
---|
2843 | ALLOCATE(lat_ful(iml)) |
---|
2844 | ALLOCATE(lon_ful(iml)) |
---|
2845 | ALLOCATE(vegmap(iml)) |
---|
2846 | ! |
---|
2847 | WRITE(numout,*) 'Reading the vegetation file' |
---|
2848 | ! |
---|
2849 | IF (is_root_prc) THEN |
---|
2850 | CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) |
---|
2851 | CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) |
---|
2852 | CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) |
---|
2853 | ! |
---|
2854 | CALL flinclo(fid) |
---|
2855 | ENDIF |
---|
2856 | |
---|
2857 | CALL bcast(lon_ful) |
---|
2858 | CALL bcast(lat_ful) |
---|
2859 | CALL bcast(vegmap) |
---|
2860 | |
---|
2861 | ! |
---|
2862 | IF (MAXVAL(vegmap) .LT. nolson) THEN |
---|
2863 | WRITE(*,*) 'WARNING -- WARNING' |
---|
2864 | WRITE(*,*) 'The vegetation map has to few vegetation types.' |
---|
2865 | WRITE(*,*) 'If you are lucky it will work but please check' |
---|
2866 | ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN |
---|
2867 | WRITE(*,*) 'More vegetation types in file than the code can' |
---|
2868 | WRITE(*,*) 'deal with.: ', MAXVAL(vegmap), nolson |
---|
2869 | STOP 'slowproc_interpol' |
---|
2870 | ENDIF |
---|
2871 | ! |
---|
2872 | ALLOCATE(lon_up(nbpt)) |
---|
2873 | ALLOCATE(lon_low(nbpt)) |
---|
2874 | ALLOCATE(lat_up(nbpt)) |
---|
2875 | ALLOCATE(lat_low(nbpt)) |
---|
2876 | ! |
---|
2877 | DO ib =1, nbpt |
---|
2878 | ! |
---|
2879 | ! We find the 4 limits of the grid-box. As we transform the resolution of the model |
---|
2880 | ! into longitudes and latitudes we do not have the problem of periodicity. |
---|
2881 | ! coslat is a help variable here ! |
---|
2882 | ! |
---|
2883 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth |
---|
2884 | ! |
---|
2885 | lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
2886 | lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
2887 | ! |
---|
2888 | coslat = pi/180. * R_Earth |
---|
2889 | ! |
---|
2890 | lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
2891 | lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
2892 | ! |
---|
2893 | ! |
---|
2894 | veget(ib,:) = 0.0 |
---|
2895 | frac_nobio (ib,:) = 0.0 |
---|
2896 | ! |
---|
2897 | ENDDO |
---|
2898 | ! |
---|
2899 | ! Get the limits of the integration domaine so that we can speed up the calculations |
---|
2900 | ! |
---|
2901 | domaine_lon_min = MINVAL(lon_low) |
---|
2902 | domaine_lon_max = MAXVAL(lon_up) |
---|
2903 | domaine_lat_min = MINVAL(lat_low) |
---|
2904 | domaine_lat_max = MAXVAL(lat_up) |
---|
2905 | ! |
---|
2906 | !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max |
---|
2907 | !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max |
---|
2908 | ! |
---|
2909 | ! Ensure that the fine grid covers the whole domain |
---|
2910 | WHERE ( lon_ful(:) .LT. domaine_lon_min ) |
---|
2911 | lon_ful(:) = lon_ful(:) + 360. |
---|
2912 | ENDWHERE |
---|
2913 | ! |
---|
2914 | WHERE ( lon_ful(:) .GT. domaine_lon_max ) |
---|
2915 | lon_ful(:) = lon_ful(:) - 360. |
---|
2916 | ENDWHERE |
---|
2917 | ! |
---|
2918 | WRITE(numout,*) 'Interpolating the vegetation map :' |
---|
2919 | WRITE(numout,'(2a40)')'0%--------------------------------------', & |
---|
2920 | & '------------------------------------100%' |
---|
2921 | ! |
---|
2922 | ilast = 1 |
---|
2923 | n_origveg(:,:) = 0 |
---|
2924 | ! |
---|
2925 | DO ip=1,iml |
---|
2926 | ! |
---|
2927 | ! Give a progress meter |
---|
2928 | ! |
---|
2929 | ! prog = ip/float(iml)*79. |
---|
2930 | ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN |
---|
2931 | ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' |
---|
2932 | ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) |
---|
2933 | ! WRITE(numout, advance="no", FMT='(a80)') meter |
---|
2934 | ! ENDIF |
---|
2935 | iprog = NINT(float(ip)/float(iml)*79.) - NINT(float(ip-1)/float(iml)*79.) |
---|
2936 | IF ( iprog .NE. 0 ) THEN |
---|
2937 | WRITE(numout,'(a1,$)') 'x' |
---|
2938 | ENDIF |
---|
2939 | ! |
---|
2940 | ! Only start looking for its place in the smaler grid if we are within the domaine |
---|
2941 | ! That should speed up things ! |
---|
2942 | ! |
---|
2943 | IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & |
---|
2944 | ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & |
---|
2945 | ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & |
---|
2946 | ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN |
---|
2947 | ! |
---|
2948 | ! look for point on GCM grid which this point on fine grid belongs to. |
---|
2949 | ! First look at the point on the model grid where we arrived just before. There is |
---|
2950 | ! a good chace that neighbouring points on the fine grid fall into the same model |
---|
2951 | ! grid box. |
---|
2952 | ! |
---|
2953 | ! |
---|
2954 | ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE |
---|
2955 | ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED. |
---|
2956 | ! |
---|
2957 | IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. & |
---|
2958 | ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & |
---|
2959 | ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. & |
---|
2960 | ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN |
---|
2961 | ! |
---|
2962 | ! We were lucky |
---|
2963 | ! |
---|
2964 | n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1 |
---|
2965 | ! |
---|
2966 | ELSE |
---|
2967 | ! |
---|
2968 | ! Otherwise, look everywhere. |
---|
2969 | ! Begin close to last grid point. |
---|
2970 | ! |
---|
2971 | found = .FALSE. |
---|
2972 | idi = 1 |
---|
2973 | ! |
---|
2974 | DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) |
---|
2975 | ! |
---|
2976 | ! forward and backward |
---|
2977 | ! |
---|
2978 | DO ii = 1,2 |
---|
2979 | ! |
---|
2980 | IF ( ii .EQ. 1 ) THEN |
---|
2981 | ib = ilast - idi |
---|
2982 | ELSE |
---|
2983 | ib = ilast + idi |
---|
2984 | ENDIF |
---|
2985 | ! |
---|
2986 | IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN |
---|
2987 | IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. & |
---|
2988 | ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & |
---|
2989 | ( lat_ful(ip) .GT. lat_low(ib) ) .AND. & |
---|
2990 | ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN |
---|
2991 | ! |
---|
2992 | n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1 |
---|
2993 | ilast = ib |
---|
2994 | found = .TRUE. |
---|
2995 | ! |
---|
2996 | ENDIF |
---|
2997 | ENDIF |
---|
2998 | ! |
---|
2999 | ENDDO |
---|
3000 | ! |
---|
3001 | idi = idi + 1 |
---|
3002 | ! |
---|
3003 | ENDDO |
---|
3004 | ! |
---|
3005 | ENDIF ! lucky/not lucky |
---|
3006 | ! |
---|
3007 | ENDIF ! in the domain |
---|
3008 | ENDDO |
---|
3009 | |
---|
3010 | ! |
---|
3011 | ! Now we know how many points of which Olson type from the fine grid fall |
---|
3012 | ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) |
---|
3013 | ! |
---|
3014 | |
---|
3015 | ! |
---|
3016 | ! determine number of points of the fine grid which fall into each box of the |
---|
3017 | ! coarse grid |
---|
3018 | ! |
---|
3019 | DO ib = 1, nbpt |
---|
3020 | n_found(ib) = SUM( n_origveg(ib,:) ) |
---|
3021 | ENDDO |
---|
3022 | |
---|
3023 | ! |
---|
3024 | ! determine fraction of Olson vegetation type in each box of the coarse grid |
---|
3025 | ! |
---|
3026 | DO vid = 1, nolson |
---|
3027 | WHERE ( n_found(:) .GT. 0 ) |
---|
3028 | frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std) |
---|
3029 | ELSEWHERE |
---|
3030 | frac_origveg(:,vid) = 0. |
---|
3031 | ENDWHERE |
---|
3032 | ENDDO |
---|
3033 | |
---|
3034 | ! |
---|
3035 | ! now finally calculate coarse vegetation map |
---|
3036 | ! Find which model vegetation corresponds to each Olson type |
---|
3037 | ! |
---|
3038 | DO vid = 1, nolson |
---|
3039 | ! |
---|
3040 | DO jv = 1, nvm |
---|
3041 | veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) |
---|
3042 | ENDDO |
---|
3043 | ! |
---|
3044 | DO jv = 1, nnobio |
---|
3045 | frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) |
---|
3046 | ENDDO |
---|
3047 | ! |
---|
3048 | ENDDO |
---|
3049 | ! |
---|
3050 | ! |
---|
3051 | WRITE(numout,*) |
---|
3052 | WRITE(numout,*) 'Interpolation Done' |
---|
3053 | ! |
---|
3054 | ! Clean up the point of the map |
---|
3055 | ! |
---|
3056 | DO ib = 1, nbpt |
---|
3057 | ! |
---|
3058 | ! Let us see if all points found something in the 5km map ! |
---|
3059 | ! |
---|
3060 | IF ( n_found(ib) .EQ. 0 ) THEN |
---|
3061 | ! |
---|
3062 | ! Now we need to handle some exceptions |
---|
3063 | ! |
---|
3064 | IF ( lalo(ib,1) .LT. -56.0) THEN |
---|
3065 | ! Antartica |
---|
3066 | frac_nobio(ib,:) = 0.0 |
---|
3067 | frac_nobio(ib,iice) = 1.0 |
---|
3068 | veget(ib,:) = 0.0 |
---|
3069 | ! |
---|
3070 | ELSE IF ( lalo(ib,1) .GT. 70.0) THEN |
---|
3071 | ! Artica |
---|
3072 | frac_nobio(ib,:) = 0.0 |
---|
3073 | frac_nobio(ib,iice) = 1.0 |
---|
3074 | veget(ib,:) = 0.0 |
---|
3075 | ! |
---|
3076 | ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN |
---|
3077 | ! Greenland |
---|
3078 | frac_nobio(ib,:) = 0.0 |
---|
3079 | frac_nobio(ib,iice) = 1.0 |
---|
3080 | veget(ib,:) = 0.0 |
---|
3081 | ! |
---|
3082 | ELSE |
---|
3083 | ! |
---|
3084 | WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box' |
---|
3085 | WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib) |
---|
3086 | WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib) |
---|
3087 | ! |
---|
3088 | WRITE(numout,*) 'Looking for nearest point on the 5 km map' |
---|
3089 | CALL slowproc_nearest (iml, lon_ful, lat_ful, & |
---|
3090 | lalo(ib,2), lalo(ib,1), inear) |
---|
3091 | WRITE(numout,*) 'Coordinates of the nearest point:', & |
---|
3092 | lon_ful(inear),lat_ful(inear) |
---|
3093 | ! |
---|
3094 | DO jv = 1, nvm |
---|
3095 | veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) |
---|
3096 | ENDDO |
---|
3097 | ! |
---|
3098 | DO jv = 1, nnobio |
---|
3099 | frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) |
---|
3100 | ENDDO |
---|
3101 | ! |
---|
3102 | ENDIF |
---|
3103 | ! |
---|
3104 | ENDIF |
---|
3105 | ! |
---|
3106 | ! |
---|
3107 | ! Limit the smalest vegetation fraction to 0.5% |
---|
3108 | ! |
---|
3109 | DO vid = 1, nvm |
---|
3110 | IF ( veget(ib,vid) .LT. min_vegfrac ) THEN |
---|
3111 | veget(ib,vid) = 0.0 |
---|
3112 | ENDIF |
---|
3113 | ENDDO |
---|
3114 | ! |
---|
3115 | sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) |
---|
3116 | frac_nobio(ib,:) = frac_nobio(ib,:)/sumf |
---|
3117 | veget(ib,:) = veget(ib,:)/sumf |
---|
3118 | ! |
---|
3119 | ! |
---|
3120 | ENDDO |
---|
3121 | ! |
---|
3122 | DEALLOCATE(lon_up) |
---|
3123 | DEALLOCATE(lon_low) |
---|
3124 | DEALLOCATE(lat_up) |
---|
3125 | DEALLOCATE(lat_low) |
---|
3126 | DEALLOCATE(lat_ful) |
---|
3127 | DEALLOCATE(lon_ful) |
---|
3128 | DEALLOCATE(vegmap) |
---|
3129 | ! |
---|
3130 | RETURN |
---|
3131 | ! |
---|
3132 | END SUBROUTINE slowproc_interpol_OLD |
---|
3133 | !! |
---|
3134 | !! Interpolate the IGBP vegetation map to the grid of the model |
---|
3135 | !! |
---|
3136 | SUBROUTINE slowproc_interpol_NEW(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio ) |
---|
3137 | ! |
---|
3138 | ! |
---|
3139 | ! |
---|
3140 | ! 0.1 INPUT |
---|
3141 | ! |
---|
3142 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
3143 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
3144 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point |
---|
3145 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
3146 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
3147 | REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid |
---|
3148 | ! |
---|
3149 | ! 0.2 OUTPUT |
---|
3150 | ! |
---|
3151 | REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions |
---|
3152 | REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... |
---|
3153 | ! |
---|
3154 | LOGICAL :: ok_interpol ! optionnal return of aggregate_vec |
---|
3155 | ! |
---|
3156 | ! 0.3 LOCAL |
---|
3157 | ! |
---|
3158 | INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes |
---|
3159 | ! |
---|
3160 | ! |
---|
3161 | CHARACTER(LEN=80) :: filename |
---|
3162 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid |
---|
3163 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap |
---|
3164 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area |
---|
3165 | INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index |
---|
3166 | REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg |
---|
3167 | REAL(r_std), DIMENSION(nbpt) :: n_found |
---|
3168 | REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg |
---|
3169 | REAL(r_std) :: vegcorr(nolson,nvm) |
---|
3170 | REAL(r_std) :: nobiocorr(nolson,nnobio) |
---|
3171 | CHARACTER(LEN=40) :: callsign |
---|
3172 | REAL(r_std) :: sumf, resol_lon, resol_lat |
---|
3173 | INTEGER(i_std) :: idi, jv, inear, nbvmax |
---|
3174 | ! |
---|
3175 | INTEGER :: ALLOC_ERR |
---|
3176 | ! |
---|
3177 | n_origveg(:,:) = zero |
---|
3178 | n_found(:) = zero |
---|
3179 | ! |
---|
3180 | CALL get_vegcorr (nolson,vegcorr,nobiocorr) |
---|
3181 | ! |
---|
3182 | !Config Key = VEGETATION_FILE |
---|
3183 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
3184 | !Config If = !IMPOSE_VEG |
---|
3185 | !Config If = !LAND_USE |
---|
3186 | !Config Def = ../surfmap/carteveg5km.nc |
---|
3187 | !Config Help = The name of the file to be opened to read the vegetation |
---|
3188 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
3189 | !Config map which is derived from the IGBP one. We assume that we have |
---|
3190 | !Config a classification in 87 types. This is Olson modified by Viovy. |
---|
3191 | ! |
---|
3192 | filename = '../surfmap/carteveg5km.nc' |
---|
3193 | CALL getin_p('VEGETATION_FILE',filename) |
---|
3194 | ! |
---|
3195 | if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
3196 | CALL bcast(iml) |
---|
3197 | CALL bcast(jml) |
---|
3198 | CALL bcast(lml) |
---|
3199 | CALL bcast(tml) |
---|
3200 | ! |
---|
3201 | ! |
---|
3202 | ALLOC_ERR=-1 |
---|
3203 | ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR) |
---|
3204 | IF (ALLOC_ERR/=0) THEN |
---|
3205 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR |
---|
3206 | STOP |
---|
3207 | ENDIF |
---|
3208 | ALLOC_ERR=-1 |
---|
3209 | ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR) |
---|
3210 | IF (ALLOC_ERR/=0) THEN |
---|
3211 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR |
---|
3212 | STOP |
---|
3213 | ENDIF |
---|
3214 | ALLOC_ERR=-1 |
---|
3215 | ALLOCATE(vegmap(iml), STAT=ALLOC_ERR) |
---|
3216 | IF (ALLOC_ERR/=0) THEN |
---|
3217 | WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR |
---|
3218 | STOP |
---|
3219 | ENDIF |
---|
3220 | ! |
---|
3221 | WRITE(numout,*) 'Reading the OLSON type vegetation file' |
---|
3222 | ! |
---|
3223 | IF (is_root_prc) THEN |
---|
3224 | CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) |
---|
3225 | CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) |
---|
3226 | CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) |
---|
3227 | ! |
---|
3228 | CALL flinclo(fid) |
---|
3229 | ENDIF |
---|
3230 | |
---|
3231 | CALL bcast(lon_ful) |
---|
3232 | CALL bcast(lat_ful) |
---|
3233 | CALL bcast(vegmap) |
---|
3234 | |
---|
3235 | ! |
---|
3236 | IF (MAXVAL(vegmap) .LT. nolson) THEN |
---|
3237 | WRITE(numout,*) 'WARNING -- WARNING' |
---|
3238 | WRITE(numout,*) 'The vegetation map has to few vegetation types.' |
---|
3239 | WRITE(numout,*) 'If you are lucky it will work but please check' |
---|
3240 | ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN |
---|
3241 | WRITE(numout,*) 'More vegetation types in file than the code can' |
---|
3242 | WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson |
---|
3243 | STOP 'slowproc_interpol' |
---|
3244 | ENDIF |
---|
3245 | ! |
---|
3246 | ! Some assumptions on the vegetation file. This information should be |
---|
3247 | ! be computed or read from the file. |
---|
3248 | ! It is the reolution in meters of the grid of the vegetation file. |
---|
3249 | ! |
---|
3250 | resol_lon = 5000. |
---|
3251 | resol_lat = 5000. |
---|
3252 | ! |
---|
3253 | ! The number of maximum vegetation map points in the GCM grid should |
---|
3254 | ! also be computed and not imposed here. |
---|
3255 | nbvmax = iml/nbpt |
---|
3256 | ! |
---|
3257 | callsign="Vegetation map" |
---|
3258 | ! |
---|
3259 | ok_interpol = .FALSE. |
---|
3260 | DO WHILE ( .NOT. ok_interpol ) |
---|
3261 | WRITE(numout,*) "Projection arrays for ",callsign," : " |
---|
3262 | WRITE(numout,*) "nbvmax = ",nbvmax |
---|
3263 | ! |
---|
3264 | ALLOC_ERR=-1 |
---|
3265 | ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
3266 | IF (ALLOC_ERR/=0) THEN |
---|
3267 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR |
---|
3268 | STOP |
---|
3269 | ENDIF |
---|
3270 | sub_index(:,:)=0 |
---|
3271 | ALLOC_ERR=-1 |
---|
3272 | ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
3273 | IF (ALLOC_ERR/=0) THEN |
---|
3274 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR |
---|
3275 | STOP |
---|
3276 | ENDIF |
---|
3277 | sub_area(:,:)=zero |
---|
3278 | ! |
---|
3279 | CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, & |
---|
3280 | & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & |
---|
3281 | & nbvmax, sub_index, sub_area, ok_interpol) |
---|
3282 | ! |
---|
3283 | IF ( .NOT. ok_interpol ) THEN |
---|
3284 | DEALLOCATE(sub_area) |
---|
3285 | DEALLOCATE(sub_index) |
---|
3286 | ! |
---|
3287 | nbvmax = nbvmax * 2 |
---|
3288 | ELSE |
---|
3289 | ! |
---|
3290 | DO ib = 1, nbpt |
---|
3291 | idi=1 |
---|
3292 | DO WHILE ( sub_area(ib,idi) > zero ) |
---|
3293 | ip = sub_index(ib,idi) |
---|
3294 | n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi) |
---|
3295 | n_found(ib) = n_found(ib) + sub_area(ib,idi) |
---|
3296 | idi = idi +1 |
---|
3297 | ENDDO |
---|
3298 | ENDDO |
---|
3299 | ! |
---|
3300 | ENDIF |
---|
3301 | ENDDO |
---|
3302 | ! |
---|
3303 | ! Now we know how many points of which Olson type from the fine grid fall |
---|
3304 | ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) |
---|
3305 | ! |
---|
3306 | ! |
---|
3307 | ! determine fraction of Olson vegetation type in each box of the coarse grid |
---|
3308 | ! |
---|
3309 | DO vid = 1, nolson |
---|
3310 | WHERE ( n_found(:) .GT. 0 ) |
---|
3311 | frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:) |
---|
3312 | ELSEWHERE |
---|
3313 | frac_origveg(:,vid) = 0. |
---|
3314 | ENDWHERE |
---|
3315 | ENDDO |
---|
3316 | ! |
---|
3317 | ! now finally calculate coarse vegetation map |
---|
3318 | ! Find which model vegetation corresponds to each Olson type |
---|
3319 | ! |
---|
3320 | veget(:,:) = zero |
---|
3321 | frac_nobio(:,:) = zero |
---|
3322 | ! |
---|
3323 | DO vid = 1, nolson |
---|
3324 | ! |
---|
3325 | DO jv = 1, nvm |
---|
3326 | veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) |
---|
3327 | ENDDO |
---|
3328 | ! |
---|
3329 | DO jv = 1, nnobio |
---|
3330 | frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) |
---|
3331 | ENDDO |
---|
3332 | ! |
---|
3333 | ENDDO |
---|
3334 | ! |
---|
3335 | WRITE(numout,*) 'slowproc_interpol : Interpolation Done' |
---|
3336 | ! |
---|
3337 | ! Clean up the point of the map |
---|
3338 | ! |
---|
3339 | DO ib = 1, nbpt |
---|
3340 | ! |
---|
3341 | ! Let us see if all points found something in the 5km map ! |
---|
3342 | ! |
---|
3343 | IF ( n_found(ib) .EQ. 0 ) THEN |
---|
3344 | ! |
---|
3345 | ! Now we need to handle some exceptions |
---|
3346 | ! |
---|
3347 | IF ( lalo(ib,1) .LT. -56.0) THEN |
---|
3348 | ! Antartica |
---|
3349 | frac_nobio(ib,:) = 0.0 |
---|
3350 | frac_nobio(ib,iice) = 1.0 |
---|
3351 | veget(ib,:) = 0.0 |
---|
3352 | ! |
---|
3353 | ELSE IF ( lalo(ib,1) .GT. 70.0) THEN |
---|
3354 | ! Artica |
---|
3355 | frac_nobio(ib,:) = 0.0 |
---|
3356 | frac_nobio(ib,iice) = 1.0 |
---|
3357 | veget(ib,:) = 0.0 |
---|
3358 | ! |
---|
3359 | ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN |
---|
3360 | ! Greenland |
---|
3361 | frac_nobio(ib,:) = 0.0 |
---|
3362 | frac_nobio(ib,iice) = 1.0 |
---|
3363 | veget(ib,:) = 0.0 |
---|
3364 | ! |
---|
3365 | ELSE |
---|
3366 | ! |
---|
3367 | WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib |
---|
3368 | WRITE(numout,*) 'Longitude range : ', lalo(ib,2) |
---|
3369 | WRITE(numout,*) 'Latitude range : ', lalo(ib,1) |
---|
3370 | ! |
---|
3371 | WRITE(numout,*) 'Looking for nearest point on the 5 km map' |
---|
3372 | CALL slowproc_nearest (iml, lon_ful, lat_ful, & |
---|
3373 | lalo(ib,2), lalo(ib,1), inear) |
---|
3374 | WRITE(numout,*) 'Coordinates of the nearest point:', & |
---|
3375 | lon_ful(inear),lat_ful(inear) |
---|
3376 | ! |
---|
3377 | DO jv = 1, nvm |
---|
3378 | veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) |
---|
3379 | ENDDO |
---|
3380 | ! |
---|
3381 | DO jv = 1, nnobio |
---|
3382 | frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) |
---|
3383 | ENDDO |
---|
3384 | ! |
---|
3385 | ENDIF |
---|
3386 | ! |
---|
3387 | ENDIF |
---|
3388 | ! |
---|
3389 | ! |
---|
3390 | ! Limit the smalest vegetation fraction to 0.5% |
---|
3391 | ! |
---|
3392 | DO vid = 1, nvm |
---|
3393 | IF ( veget(ib,vid) .LT. min_vegfrac ) THEN |
---|
3394 | veget(ib,vid) = 0.0 |
---|
3395 | ENDIF |
---|
3396 | ENDDO |
---|
3397 | ! |
---|
3398 | sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) |
---|
3399 | frac_nobio(ib,:) = frac_nobio(ib,:)/sumf |
---|
3400 | veget(ib,:) = veget(ib,:)/sumf |
---|
3401 | ! |
---|
3402 | ! |
---|
3403 | ENDDO |
---|
3404 | ! |
---|
3405 | DEALLOCATE(vegmap) |
---|
3406 | DEALLOCATE(lat_ful, lon_ful) |
---|
3407 | DEALLOCATE(sub_index) |
---|
3408 | DEALLOCATE(sub_area) |
---|
3409 | |
---|
3410 | ! |
---|
3411 | RETURN |
---|
3412 | ! |
---|
3413 | END SUBROUTINE slowproc_interpol_NEW |
---|
3414 | |
---|
3415 | !! |
---|
3416 | !! Interpolate the IGBP vegetation map to the grid of the model |
---|
3417 | !MM TAG 1.6 model ! |
---|
3418 | !! |
---|
3419 | SUBROUTINE slowproc_interpol_OLD_g(nbpt, lalo, neighbours, resolution, veget, frac_nobio ) |
---|
3420 | ! |
---|
3421 | ! |
---|
3422 | ! |
---|
3423 | ! 0.1 INPUT |
---|
3424 | ! |
---|
3425 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
3426 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
3427 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point |
---|
3428 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
3429 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
3430 | ! |
---|
3431 | ! 0.2 OUTPUT |
---|
3432 | ! |
---|
3433 | REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions |
---|
3434 | REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... |
---|
3435 | ! |
---|
3436 | ! 0.3 LOCAL |
---|
3437 | ! |
---|
3438 | REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8 |
---|
3439 | INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes |
---|
3440 | ! |
---|
3441 | ! |
---|
3442 | CHARACTER(LEN=80) :: filename |
---|
3443 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid |
---|
3444 | REAL(r_std) :: lev(1), date, dt, coslat, pi |
---|
3445 | INTEGER(i_std) :: itau(1) |
---|
3446 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap |
---|
3447 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low |
---|
3448 | INTEGER, DIMENSION(nbpt,nolson) :: n_origveg |
---|
3449 | INTEGER, DIMENSION(nbpt) :: n_found |
---|
3450 | REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg |
---|
3451 | REAL(r_std) :: vegcorr(nolson,nvm) |
---|
3452 | REAL(r_std) :: nobiocorr(nolson,nnobio) |
---|
3453 | CHARACTER(LEN=80) :: meter |
---|
3454 | REAL(r_std) :: prog, sumf |
---|
3455 | LOGICAL :: found |
---|
3456 | INTEGER :: idi, ilast, ii, jv, inear, iprog |
---|
3457 | REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max |
---|
3458 | ! |
---|
3459 | pi = 4. * ATAN(1.) |
---|
3460 | ! |
---|
3461 | CALL get_vegcorr (nolson,vegcorr,nobiocorr) |
---|
3462 | ! |
---|
3463 | !Config Key = VEGETATION_FILE |
---|
3464 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
3465 | !Config If = !IMPOSE_VEG |
---|
3466 | !Config Def = ../surfmap/carteveg5km.nc |
---|
3467 | !Config Help = The name of the file to be opened to read the vegetation |
---|
3468 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
3469 | !Config map which is derived from the IGBP one. We assume that we have |
---|
3470 | !Config a classification in 87 types. This is Olson modified by Viovy. |
---|
3471 | ! |
---|
3472 | filename = '../surfmap/carteveg5km.nc' |
---|
3473 | CALL getin('VEGETATION_FILE',filename) |
---|
3474 | ! |
---|
3475 | CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
3476 | ! |
---|
3477 | ! |
---|
3478 | ALLOCATE(lat_ful(iml)) |
---|
3479 | ALLOCATE(lon_ful(iml)) |
---|
3480 | ALLOCATE(vegmap(iml)) |
---|
3481 | ! |
---|
3482 | WRITE(numout,*) 'Reading the vegetation file' |
---|
3483 | ! |
---|
3484 | CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) |
---|
3485 | CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) |
---|
3486 | CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) |
---|
3487 | ! |
---|
3488 | CALL flinclo(fid) |
---|
3489 | |
---|
3490 | ! |
---|
3491 | IF (MAXVAL(vegmap) .LT. nolson) THEN |
---|
3492 | WRITE(*,*) 'WARNING -- WARNING' |
---|
3493 | WRITE(*,*) 'The vegetation map has to few vegetation types.' |
---|
3494 | WRITE(*,*) 'If you are lucky it will work but please check' |
---|
3495 | ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN |
---|
3496 | WRITE(*,*) 'More vegetation types in file than the code can' |
---|
3497 | WRITE(*,*) 'deal with.: ', MAXVAL(vegmap), nolson |
---|
3498 | STOP 'slowproc_interpol' |
---|
3499 | ENDIF |
---|
3500 | ! |
---|
3501 | ALLOCATE(lon_up(nbpt)) |
---|
3502 | ALLOCATE(lon_low(nbpt)) |
---|
3503 | ALLOCATE(lat_up(nbpt)) |
---|
3504 | ALLOCATE(lat_low(nbpt)) |
---|
3505 | ! |
---|
3506 | DO ib =1, nbpt |
---|
3507 | ! |
---|
3508 | ! We find the 4 limits of the grid-box. As we transform the resolution of the model |
---|
3509 | ! into longitudes and latitudes we do not have the problem of periodicity. |
---|
3510 | ! coslat is a help variable here ! |
---|
3511 | ! |
---|
3512 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth |
---|
3513 | ! |
---|
3514 | lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
3515 | lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
3516 | ! |
---|
3517 | coslat = pi/180. * R_Earth |
---|
3518 | ! |
---|
3519 | lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
3520 | lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
3521 | ! |
---|
3522 | ! |
---|
3523 | veget(ib,:) = 0.0 |
---|
3524 | frac_nobio (ib,:) = 0.0 |
---|
3525 | ! |
---|
3526 | ENDDO |
---|
3527 | ! |
---|
3528 | ! Get the limits of the integration domaine so that we can speed up the calculations |
---|
3529 | ! |
---|
3530 | domaine_lon_min = MINVAL(lon_low) |
---|
3531 | domaine_lon_max = MAXVAL(lon_up) |
---|
3532 | domaine_lat_min = MINVAL(lat_low) |
---|
3533 | domaine_lat_max = MAXVAL(lat_up) |
---|
3534 | ! |
---|
3535 | !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max |
---|
3536 | !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max |
---|
3537 | ! |
---|
3538 | ! Ensure that the fine grid covers the whole domain |
---|
3539 | WHERE ( lon_ful(:) .LT. domaine_lon_min ) |
---|
3540 | lon_ful(:) = lon_ful(:) + 360. |
---|
3541 | ENDWHERE |
---|
3542 | ! |
---|
3543 | WHERE ( lon_ful(:) .GT. domaine_lon_max ) |
---|
3544 | lon_ful(:) = lon_ful(:) - 360. |
---|
3545 | ENDWHERE |
---|
3546 | ! |
---|
3547 | WRITE(numout,*) 'Interpolating the vegetation map :' |
---|
3548 | WRITE(numout,'(2a40)')'0%--------------------------------------', & |
---|
3549 | & '------------------------------------100%' |
---|
3550 | ! |
---|
3551 | ilast = 1 |
---|
3552 | n_origveg(:,:) = 0 |
---|
3553 | ! |
---|
3554 | DO ip=1,iml |
---|
3555 | ! |
---|
3556 | ! Give a progress meter |
---|
3557 | ! |
---|
3558 | ! prog = ip/float(iml)*79. |
---|
3559 | ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN |
---|
3560 | ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' |
---|
3561 | ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) |
---|
3562 | ! WRITE(numout, advance="no", FMT='(a80)') meter |
---|
3563 | ! ENDIF |
---|
3564 | iprog = NINT(float(ip)/float(iml)*79.) - NINT(float(ip-1)/float(iml)*79.) |
---|
3565 | IF ( iprog .NE. 0 ) THEN |
---|
3566 | WRITE(numout,'(a1,$)') 'x' |
---|
3567 | ENDIF |
---|
3568 | ! |
---|
3569 | ! Only start looking for its place in the smaler grid if we are within the domaine |
---|
3570 | ! That should speed up things ! |
---|
3571 | ! |
---|
3572 | IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & |
---|
3573 | ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & |
---|
3574 | ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & |
---|
3575 | ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN |
---|
3576 | ! |
---|
3577 | ! look for point on GCM grid which this point on fine grid belongs to. |
---|
3578 | ! First look at the point on the model grid where we arrived just before. There is |
---|
3579 | ! a good chace that neighbouring points on the fine grid fall into the same model |
---|
3580 | ! grid box. |
---|
3581 | ! |
---|
3582 | ! |
---|
3583 | ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE |
---|
3584 | ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED. |
---|
3585 | ! |
---|
3586 | IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. & |
---|
3587 | ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & |
---|
3588 | ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. & |
---|
3589 | ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN |
---|
3590 | ! |
---|
3591 | ! We were lucky |
---|
3592 | ! |
---|
3593 | n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1 |
---|
3594 | ! |
---|
3595 | ELSE |
---|
3596 | ! |
---|
3597 | ! Otherwise, look everywhere. |
---|
3598 | ! Begin close to last grid point. |
---|
3599 | ! |
---|
3600 | found = .FALSE. |
---|
3601 | idi = 1 |
---|
3602 | ! |
---|
3603 | DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) |
---|
3604 | ! |
---|
3605 | ! forward and backward |
---|
3606 | ! |
---|
3607 | DO ii = 1,2 |
---|
3608 | ! |
---|
3609 | IF ( ii .EQ. 1 ) THEN |
---|
3610 | ib = ilast - idi |
---|
3611 | ELSE |
---|
3612 | ib = ilast + idi |
---|
3613 | ENDIF |
---|
3614 | ! |
---|
3615 | IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN |
---|
3616 | IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. & |
---|
3617 | ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & |
---|
3618 | ( lat_ful(ip) .GT. lat_low(ib) ) .AND. & |
---|
3619 | ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN |
---|
3620 | ! |
---|
3621 | n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1 |
---|
3622 | ilast = ib |
---|
3623 | found = .TRUE. |
---|
3624 | ! |
---|
3625 | ENDIF |
---|
3626 | ENDIF |
---|
3627 | ! |
---|
3628 | ENDDO |
---|
3629 | ! |
---|
3630 | idi = idi + 1 |
---|
3631 | ! |
---|
3632 | ENDDO |
---|
3633 | ! |
---|
3634 | ENDIF ! lucky/not lucky |
---|
3635 | ! |
---|
3636 | ENDIF ! in the domain |
---|
3637 | ENDDO |
---|
3638 | |
---|
3639 | ! |
---|
3640 | ! Now we know how many points of which Olson type from the fine grid fall |
---|
3641 | ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) |
---|
3642 | ! |
---|
3643 | |
---|
3644 | ! |
---|
3645 | ! determine number of points of the fine grid which fall into each box of the |
---|
3646 | ! coarse grid |
---|
3647 | ! |
---|
3648 | DO ib = 1, nbpt |
---|
3649 | n_found(ib) = SUM( n_origveg(ib,:) ) |
---|
3650 | ENDDO |
---|
3651 | |
---|
3652 | ! |
---|
3653 | ! determine fraction of Olson vegetation type in each box of the coarse grid |
---|
3654 | ! |
---|
3655 | DO vid = 1, nolson |
---|
3656 | WHERE ( n_found(:) .GT. 0 ) |
---|
3657 | frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std) |
---|
3658 | ELSEWHERE |
---|
3659 | frac_origveg(:,vid) = 0. |
---|
3660 | ENDWHERE |
---|
3661 | ENDDO |
---|
3662 | |
---|
3663 | ! |
---|
3664 | ! now finally calculate coarse vegetation map |
---|
3665 | ! Find which model vegetation corresponds to each Olson type |
---|
3666 | ! |
---|
3667 | DO vid = 1, nolson |
---|
3668 | ! |
---|
3669 | DO jv = 1, nvm |
---|
3670 | veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) |
---|
3671 | ENDDO |
---|
3672 | ! |
---|
3673 | DO jv = 1, nnobio |
---|
3674 | frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) |
---|
3675 | ENDDO |
---|
3676 | ! |
---|
3677 | ENDDO |
---|
3678 | ! |
---|
3679 | ! |
---|
3680 | WRITE(numout,*) |
---|
3681 | WRITE(numout,*) 'Interpolation Done' |
---|
3682 | ! |
---|
3683 | ! Clean up the point of the map |
---|
3684 | ! |
---|
3685 | DO ib = 1, nbpt |
---|
3686 | ! |
---|
3687 | ! Let us see if all points found something in the 5km map ! |
---|
3688 | ! |
---|
3689 | IF ( n_found(ib) .EQ. 0 ) THEN |
---|
3690 | ! |
---|
3691 | ! Now we need to handle some exceptions |
---|
3692 | ! |
---|
3693 | IF ( lalo(ib,1) .LT. -56.0) THEN |
---|
3694 | ! Antartica |
---|
3695 | frac_nobio(ib,:) = 0.0 |
---|
3696 | frac_nobio(ib,iice) = 1.0 |
---|
3697 | veget(ib,:) = 0.0 |
---|
3698 | ! |
---|
3699 | ELSE IF ( lalo(ib,1) .GT. 70.0) THEN |
---|
3700 | ! Artica |
---|
3701 | frac_nobio(ib,:) = 0.0 |
---|
3702 | frac_nobio(ib,iice) = 1.0 |
---|
3703 | veget(ib,:) = 0.0 |
---|
3704 | ! |
---|
3705 | ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN |
---|
3706 | ! Greenland |
---|
3707 | frac_nobio(ib,:) = 0.0 |
---|
3708 | frac_nobio(ib,iice) = 1.0 |
---|
3709 | veget(ib,:) = 0.0 |
---|
3710 | ! |
---|
3711 | ELSE |
---|
3712 | ! |
---|
3713 | WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box' |
---|
3714 | WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib) |
---|
3715 | WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib) |
---|
3716 | ! |
---|
3717 | WRITE(numout,*) 'Looking for nearest point on the 5 km map' |
---|
3718 | CALL slowproc_nearest (iml, lon_ful, lat_ful, & |
---|
3719 | lalo(ib,2), lalo(ib,1), inear) |
---|
3720 | WRITE(numout,*) 'Coordinates of the nearest point:', & |
---|
3721 | lon_ful(inear),lat_ful(inear) |
---|
3722 | ! |
---|
3723 | DO jv = 1, nvm |
---|
3724 | veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) |
---|
3725 | ENDDO |
---|
3726 | ! |
---|
3727 | DO jv = 1, nnobio |
---|
3728 | frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) |
---|
3729 | ENDDO |
---|
3730 | ! |
---|
3731 | ENDIF |
---|
3732 | ! |
---|
3733 | ENDIF |
---|
3734 | ! |
---|
3735 | ! |
---|
3736 | ! Limit the smalest vegetation fraction to 0.5% |
---|
3737 | ! |
---|
3738 | DO vid = 1, nvm |
---|
3739 | IF ( veget(ib,vid) .LT. min_vegfrac ) THEN |
---|
3740 | veget(ib,vid) = 0.0 |
---|
3741 | ENDIF |
---|
3742 | ENDDO |
---|
3743 | ! |
---|
3744 | sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) |
---|
3745 | frac_nobio(ib,:) = frac_nobio(ib,:)/sumf |
---|
3746 | veget(ib,:) = veget(ib,:)/sumf |
---|
3747 | ! |
---|
3748 | ! |
---|
3749 | ENDDO |
---|
3750 | ! |
---|
3751 | DEALLOCATE(lon_up) |
---|
3752 | DEALLOCATE(lon_low) |
---|
3753 | DEALLOCATE(lat_up) |
---|
3754 | DEALLOCATE(lat_low) |
---|
3755 | DEALLOCATE(lat_ful) |
---|
3756 | DEALLOCATE(lon_ful) |
---|
3757 | DEALLOCATE(vegmap) |
---|
3758 | ! |
---|
3759 | RETURN |
---|
3760 | ! |
---|
3761 | END SUBROUTINE slowproc_interpol_OLD_g |
---|
3762 | !! |
---|
3763 | !! Interpolate the IGBP vegetation map to the grid of the model |
---|
3764 | !! |
---|
3765 | SUBROUTINE slowproc_interpol_NEW_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio ) |
---|
3766 | ! |
---|
3767 | ! |
---|
3768 | ! |
---|
3769 | ! 0.1 INPUT |
---|
3770 | ! |
---|
3771 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
3772 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
3773 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point |
---|
3774 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
3775 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
3776 | REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid |
---|
3777 | ! |
---|
3778 | ! 0.2 OUTPUT |
---|
3779 | ! |
---|
3780 | REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions |
---|
3781 | REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... |
---|
3782 | ! |
---|
3783 | LOGICAL :: ok_interpol ! optionnal return of aggregate_vec |
---|
3784 | ! |
---|
3785 | ! 0.3 LOCAL |
---|
3786 | ! |
---|
3787 | INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes |
---|
3788 | ! |
---|
3789 | ! |
---|
3790 | CHARACTER(LEN=80) :: filename |
---|
3791 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid |
---|
3792 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap |
---|
3793 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area |
---|
3794 | INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index |
---|
3795 | REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg |
---|
3796 | REAL(r_std), DIMENSION(nbpt) :: n_found |
---|
3797 | REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg |
---|
3798 | REAL(r_std) :: vegcorr(nolson,nvm) |
---|
3799 | REAL(r_std) :: nobiocorr(nolson,nnobio) |
---|
3800 | CHARACTER(LEN=40) :: callsign |
---|
3801 | REAL(r_std) :: sumf, resol_lon, resol_lat |
---|
3802 | INTEGER(i_std) :: idi, jv, inear, nbvmax |
---|
3803 | ! |
---|
3804 | INTEGER :: ALLOC_ERR |
---|
3805 | ! |
---|
3806 | n_origveg(:,:) = zero |
---|
3807 | n_found(:) = zero |
---|
3808 | ! |
---|
3809 | CALL get_vegcorr (nolson,vegcorr,nobiocorr) |
---|
3810 | ! |
---|
3811 | !Config Key = VEGETATION_FILE |
---|
3812 | !Config Desc = Name of file from which the vegetation map is to be read |
---|
3813 | !Config If = !IMPOSE_VEG |
---|
3814 | !Config If = !LAND_USE |
---|
3815 | !Config Def = ../surfmap/carteveg5km.nc |
---|
3816 | !Config Help = The name of the file to be opened to read the vegetation |
---|
3817 | !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km |
---|
3818 | !Config map which is derived from the IGBP one. We assume that we have |
---|
3819 | !Config a classification in 87 types. This is Olson modified by Viovy. |
---|
3820 | ! |
---|
3821 | filename = '../surfmap/carteveg5km.nc' |
---|
3822 | CALL getin('VEGETATION_FILE',filename) |
---|
3823 | ! |
---|
3824 | CALL flininfo(filename, iml, jml, lml, tml, fid) |
---|
3825 | ! |
---|
3826 | ! |
---|
3827 | ALLOC_ERR=-1 |
---|
3828 | ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR) |
---|
3829 | IF (ALLOC_ERR/=0) THEN |
---|
3830 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR |
---|
3831 | STOP |
---|
3832 | ENDIF |
---|
3833 | ALLOC_ERR=-1 |
---|
3834 | ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR) |
---|
3835 | IF (ALLOC_ERR/=0) THEN |
---|
3836 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR |
---|
3837 | STOP |
---|
3838 | ENDIF |
---|
3839 | ALLOC_ERR=-1 |
---|
3840 | ALLOCATE(vegmap(iml), STAT=ALLOC_ERR) |
---|
3841 | IF (ALLOC_ERR/=0) THEN |
---|
3842 | WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR |
---|
3843 | STOP |
---|
3844 | ENDIF |
---|
3845 | ! |
---|
3846 | WRITE(numout,*) 'Reading the OLSON type vegetation file' |
---|
3847 | ! |
---|
3848 | CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) |
---|
3849 | CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) |
---|
3850 | CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) |
---|
3851 | ! |
---|
3852 | CALL flinclo(fid) |
---|
3853 | ! |
---|
3854 | IF (MAXVAL(vegmap) .LT. nolson) THEN |
---|
3855 | WRITE(numout,*) 'WARNING -- WARNING' |
---|
3856 | WRITE(numout,*) 'The vegetation map has to few vegetation types.' |
---|
3857 | WRITE(numout,*) 'If you are lucky it will work but please check' |
---|
3858 | ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN |
---|
3859 | WRITE(numout,*) 'More vegetation types in file than the code can' |
---|
3860 | WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson |
---|
3861 | STOP 'slowproc_interpol' |
---|
3862 | ENDIF |
---|
3863 | ! |
---|
3864 | ! Some assumptions on the vegetation file. This information should be |
---|
3865 | ! be computed or read from the file. |
---|
3866 | ! It is the reolution in meters of the grid of the vegetation file. |
---|
3867 | ! |
---|
3868 | resol_lon = 5000. |
---|
3869 | resol_lat = 5000. |
---|
3870 | ! |
---|
3871 | ! The number of maximum vegetation map points in the GCM grid should |
---|
3872 | ! also be computed and not imposed here. |
---|
3873 | nbvmax = iml/nbpt |
---|
3874 | ! |
---|
3875 | callsign="Vegetation map" |
---|
3876 | ! |
---|
3877 | ok_interpol = .FALSE. |
---|
3878 | DO WHILE ( .NOT. ok_interpol ) |
---|
3879 | WRITE(numout,*) "Projection arrays for ",callsign," : " |
---|
3880 | WRITE(numout,*) "nbvmax = ",nbvmax |
---|
3881 | ! |
---|
3882 | ALLOC_ERR=-1 |
---|
3883 | ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
3884 | IF (ALLOC_ERR/=0) THEN |
---|
3885 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR |
---|
3886 | STOP |
---|
3887 | ENDIF |
---|
3888 | sub_index(:,:)=0 |
---|
3889 | ALLOC_ERR=-1 |
---|
3890 | ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) |
---|
3891 | IF (ALLOC_ERR/=0) THEN |
---|
3892 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR |
---|
3893 | STOP |
---|
3894 | ENDIF |
---|
3895 | sub_area(:,:)=zero |
---|
3896 | ! |
---|
3897 | CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, & |
---|
3898 | & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & |
---|
3899 | & nbvmax, sub_index, sub_area, ok_interpol) |
---|
3900 | ! |
---|
3901 | IF ( .NOT. ok_interpol ) THEN |
---|
3902 | DEALLOCATE(sub_area) |
---|
3903 | DEALLOCATE(sub_index) |
---|
3904 | ! |
---|
3905 | nbvmax = nbvmax * 2 |
---|
3906 | ELSE |
---|
3907 | ! |
---|
3908 | DO ib = 1, nbpt |
---|
3909 | idi=1 |
---|
3910 | DO WHILE ( sub_area(ib,idi) > zero ) |
---|
3911 | ip = sub_index(ib,idi) |
---|
3912 | n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi) |
---|
3913 | n_found(ib) = n_found(ib) + sub_area(ib,idi) |
---|
3914 | idi = idi +1 |
---|
3915 | ENDDO |
---|
3916 | ENDDO |
---|
3917 | ! |
---|
3918 | ENDIF |
---|
3919 | ENDDO |
---|
3920 | ! |
---|
3921 | ! Now we know how many points of which Olson type from the fine grid fall |
---|
3922 | ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) |
---|
3923 | ! |
---|
3924 | ! |
---|
3925 | ! determine fraction of Olson vegetation type in each box of the coarse grid |
---|
3926 | ! |
---|
3927 | DO vid = 1, nolson |
---|
3928 | WHERE ( n_found(:) .GT. 0 ) |
---|
3929 | frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:) |
---|
3930 | ELSEWHERE |
---|
3931 | frac_origveg(:,vid) = 0. |
---|
3932 | ENDWHERE |
---|
3933 | ENDDO |
---|
3934 | ! |
---|
3935 | ! now finally calculate coarse vegetation map |
---|
3936 | ! Find which model vegetation corresponds to each Olson type |
---|
3937 | ! |
---|
3938 | veget(:,:) = zero |
---|
3939 | frac_nobio(:,:) = zero |
---|
3940 | ! |
---|
3941 | DO vid = 1, nolson |
---|
3942 | ! |
---|
3943 | DO jv = 1, nvm |
---|
3944 | veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) |
---|
3945 | ENDDO |
---|
3946 | ! |
---|
3947 | DO jv = 1, nnobio |
---|
3948 | frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) |
---|
3949 | ENDDO |
---|
3950 | ! |
---|
3951 | ENDDO |
---|
3952 | ! |
---|
3953 | WRITE(numout,*) 'slowproc_interpol : Interpolation Done' |
---|
3954 | ! |
---|
3955 | ! Clean up the point of the map |
---|
3956 | ! |
---|
3957 | DO ib = 1, nbpt |
---|
3958 | ! |
---|
3959 | ! Let us see if all points found something in the 5km map ! |
---|
3960 | ! |
---|
3961 | IF ( n_found(ib) .EQ. 0 ) THEN |
---|
3962 | ! |
---|
3963 | ! Now we need to handle some exceptions |
---|
3964 | ! |
---|
3965 | IF ( lalo(ib,1) .LT. -56.0) THEN |
---|
3966 | ! Antartica |
---|
3967 | frac_nobio(ib,:) = 0.0 |
---|
3968 | frac_nobio(ib,iice) = 1.0 |
---|
3969 | veget(ib,:) = 0.0 |
---|
3970 | ! |
---|
3971 | ELSE IF ( lalo(ib,1) .GT. 70.0) THEN |
---|
3972 | ! Artica |
---|
3973 | frac_nobio(ib,:) = 0.0 |
---|
3974 | frac_nobio(ib,iice) = 1.0 |
---|
3975 | veget(ib,:) = 0.0 |
---|
3976 | ! |
---|
3977 | ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN |
---|
3978 | ! Greenland |
---|
3979 | frac_nobio(ib,:) = 0.0 |
---|
3980 | frac_nobio(ib,iice) = 1.0 |
---|
3981 | veget(ib,:) = 0.0 |
---|
3982 | ! |
---|
3983 | ELSE |
---|
3984 | ! |
---|
3985 | WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib |
---|
3986 | WRITE(numout,*) 'Longitude range : ', lalo(ib,2) |
---|
3987 | WRITE(numout,*) 'Latitude range : ', lalo(ib,1) |
---|
3988 | ! |
---|
3989 | WRITE(numout,*) 'Looking for nearest point on the 5 km map' |
---|
3990 | CALL slowproc_nearest (iml, lon_ful, lat_ful, & |
---|
3991 | lalo(ib,2), lalo(ib,1), inear) |
---|
3992 | WRITE(numout,*) 'Coordinates of the nearest point:', & |
---|
3993 | lon_ful(inear),lat_ful(inear) |
---|
3994 | ! |
---|
3995 | DO jv = 1, nvm |
---|
3996 | veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) |
---|
3997 | ENDDO |
---|
3998 | ! |
---|
3999 | DO jv = 1, nnobio |
---|
4000 | frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) |
---|
4001 | ENDDO |
---|
4002 | ! |
---|
4003 | ENDIF |
---|
4004 | ! |
---|
4005 | ENDIF |
---|
4006 | ! |
---|
4007 | ! |
---|
4008 | ! Limit the smalest vegetation fraction to 0.5% |
---|
4009 | ! |
---|
4010 | DO vid = 1, nvm |
---|
4011 | IF ( veget(ib,vid) .LT. min_vegfrac ) THEN |
---|
4012 | veget(ib,vid) = 0.0 |
---|
4013 | ENDIF |
---|
4014 | ENDDO |
---|
4015 | ! |
---|
4016 | sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) |
---|
4017 | frac_nobio(ib,:) = frac_nobio(ib,:)/sumf |
---|
4018 | veget(ib,:) = veget(ib,:)/sumf |
---|
4019 | ! |
---|
4020 | ! |
---|
4021 | ENDDO |
---|
4022 | ! |
---|
4023 | DEALLOCATE(vegmap) |
---|
4024 | DEALLOCATE(lat_ful, lon_ful) |
---|
4025 | DEALLOCATE(sub_index) |
---|
4026 | DEALLOCATE(sub_area) |
---|
4027 | |
---|
4028 | ! |
---|
4029 | RETURN |
---|
4030 | ! |
---|
4031 | END SUBROUTINE slowproc_interpol_NEW_g |
---|
4032 | |
---|
4033 | |
---|
4034 | !! |
---|
4035 | !! looks for nearest grid point on the fine map |
---|
4036 | !! |
---|
4037 | SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear) |
---|
4038 | |
---|
4039 | INTEGER(i_std), INTENT(in) :: iml |
---|
4040 | REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 |
---|
4041 | REAL(r_std), INTENT(in) :: lonmod, latmod |
---|
4042 | |
---|
4043 | INTEGER(i_std), INTENT(out) :: inear |
---|
4044 | |
---|
4045 | REAL(r_std) :: pi |
---|
4046 | REAL(r_std) :: pa, p |
---|
4047 | REAL(r_std) :: coscolat, sincolat |
---|
4048 | REAL(r_std) :: cospa, sinpa |
---|
4049 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang |
---|
4050 | INTEGER(i_std) :: i |
---|
4051 | INTEGER(i_std), DIMENSION(1) :: ineartab |
---|
4052 | INTEGER :: ALLOC_ERR |
---|
4053 | |
---|
4054 | ALLOC_ERR=-1 |
---|
4055 | ALLOCATE(cosang(iml), STAT=ALLOC_ERR) |
---|
4056 | IF (ALLOC_ERR/=0) THEN |
---|
4057 | WRITE(numout,*) "ERROR IN ALLOCATION of cosang : ",ALLOC_ERR |
---|
4058 | STOP |
---|
4059 | ENDIF |
---|
4060 | |
---|
4061 | pi = 4.0 * ATAN(1.0) |
---|
4062 | |
---|
4063 | pa = pi/2.0 - latmod*pi/180.0 ! dist. entre pole n et point a |
---|
4064 | cospa = COS(pa) |
---|
4065 | sinpa = SIN(pa) |
---|
4066 | |
---|
4067 | DO i = 1, iml |
---|
4068 | |
---|
4069 | sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) |
---|
4070 | coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) |
---|
4071 | |
---|
4072 | p = (lonmod-lon5(i))*pi/180.0 ! angle entre a et b (leurs meridiens) |
---|
4073 | |
---|
4074 | ! dist(i) = ACOS( cospa*coscolat + sinpa*sincolat*COS(p)) |
---|
4075 | cosang(i) = cospa*coscolat + sinpa*sincolat*COS(p) |
---|
4076 | |
---|
4077 | ENDDO |
---|
4078 | |
---|
4079 | ineartab = MAXLOC( cosang(:) ) |
---|
4080 | inear = ineartab(1) |
---|
4081 | |
---|
4082 | DEALLOCATE(cosang) |
---|
4083 | END SUBROUTINE slowproc_nearest |
---|
4084 | |
---|
4085 | !! |
---|
4086 | !! Interpolate the Zobler soil type map |
---|
4087 | !! |
---|
4088 | SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soiltype, clayfraction) |
---|
4089 | ! |
---|
4090 | ! |
---|
4091 | ! This subroutine should read the Zobler map and interpolate to the model grid. The method |
---|
4092 | ! is to get fraction of the three main soiltypes for each grid box. |
---|
4093 | ! The soil fraction are going to be put into the array soiltype in the following order : |
---|
4094 | ! coarse, medium and fine. |
---|
4095 | ! |
---|
4096 | ! |
---|
4097 | ! 0.1 INPUT |
---|
4098 | ! |
---|
4099 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
4100 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
4101 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point |
---|
4102 | ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) |
---|
4103 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
4104 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
4105 | ! |
---|
4106 | ! 0.2 OUTPUT |
---|
4107 | ! |
---|
4108 | REAL(r_std), INTENT(out) :: soiltype(nbpt, nstm) ! Soil type map to be created from the Zobler map |
---|
4109 | REAL(r_std), INTENT(out) :: clayfraction(nbpt) ! The fraction of clay as used by STOMATE |
---|
4110 | ! |
---|
4111 | ! |
---|
4112 | ! 0.3 LOCAL |
---|
4113 | ! |
---|
4114 | INTEGER(i_std) :: nbvmax |
---|
4115 | ! |
---|
4116 | CHARACTER(LEN=80) :: filename |
---|
4117 | INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, nbexp |
---|
4118 | REAL(r_std) :: lev(1), date, dt |
---|
4119 | INTEGER(i_std) :: itau(1) |
---|
4120 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, soiltext |
---|
4121 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask |
---|
4122 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area |
---|
4123 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index |
---|
4124 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt |
---|
4125 | REAL(r_std) :: sgn |
---|
4126 | CHARACTER(LEN=30) :: callsign |
---|
4127 | ! |
---|
4128 | ! Number of texture classes in Zobler |
---|
4129 | ! |
---|
4130 | INTEGER(i_std), PARAMETER :: classnb = 7 |
---|
4131 | REAL(r_std) :: textfrac_table(classnb, nstm) |
---|
4132 | ! |
---|
4133 | LOGICAL :: ok_interpol ! optionnal return of aggregate_2d |
---|
4134 | ! |
---|
4135 | INTEGER :: ALLOC_ERR |
---|
4136 | ! |
---|
4137 | ! |
---|
4138 | CALL get_soilcorr (classnb, textfrac_table) |
---|
4139 | ! |
---|
4140 | ! Needs to be a configurable variable |
---|
4141 | ! |
---|
4142 | ! |
---|
4143 | !Config Key = SOILTYPE_FILE |
---|
4144 | !Config Desc = Name of file from which soil types are read |
---|
4145 | !Config Def = ../surfmap/soils_param.nc |
---|
4146 | !Config If = !IMPOSE_VEG |
---|
4147 | !Config Help = The name of the file to be opened to read the soil types. |
---|
4148 | !Config The data from this file is then interpolated to the grid of |
---|
4149 | !Config of the model. The aim is to get fractions for sand loam and |
---|
4150 | !Config clay in each grid box. This information is used for soil hydrology |
---|
4151 | !Config and respiration. |
---|
4152 | ! |
---|
4153 | filename = '../surfmap/soils_param.nc' |
---|
4154 | CALL getin_p('SOILTYPE_FILE',filename) |
---|
4155 | ! |
---|
4156 | IF (is_root_prc) THEN |
---|
4157 | CALL flininfo(filename,iml, jml, lml, tml, fid) |
---|
4158 | CALL flinclo(fid) |
---|
4159 | ENDIF |
---|
4160 | CALL bcast(iml) |
---|
4161 | CALL bcast(jml) |
---|
4162 | CALL bcast(lml) |
---|
4163 | CALL bcast(tml) |
---|
4164 | ! |
---|
4165 | ! soils_param.nc file is 1° soit texture file. |
---|
4166 | ! |
---|
4167 | ALLOC_ERR=-1 |
---|
4168 | ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR) |
---|
4169 | IF (ALLOC_ERR/=0) THEN |
---|
4170 | WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR |
---|
4171 | STOP |
---|
4172 | ENDIF |
---|
4173 | ALLOC_ERR=-1 |
---|
4174 | ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR) |
---|
4175 | IF (ALLOC_ERR/=0) THEN |
---|
4176 | WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR |
---|
4177 | STOP |
---|
4178 | ENDIF |
---|
4179 | ALLOC_ERR=-1 |
---|
4180 | ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) |
---|
4181 | IF (ALLOC_ERR/=0) THEN |
---|
4182 | WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR |
---|
4183 | STOP |
---|
4184 | ENDIF |
---|
4185 | ALLOC_ERR=-1 |
---|
4186 | ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR) |
---|
4187 | IF (ALLOC_ERR/=0) THEN |
---|
4188 | WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR |
---|
4189 | STOP |
---|
4190 | ENDIF |
---|
4191 | ! |
---|
4192 | IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid) |
---|
4193 | CALL bcast(lon_rel) |
---|
4194 | CALL bcast(lat_rel) |
---|
4195 | CALL bcast(itau) |
---|
4196 | CALL bcast(date) |
---|
4197 | CALL bcast(dt) |
---|
4198 | |
---|
4199 | ! |
---|
4200 | IF (is_root_prc) CALL flinget(fid, 'soiltext', iml, jml, lml, tml, 1, 1, soiltext) |
---|
4201 | CALL bcast(soiltext) |
---|
4202 | ! |
---|
4203 | IF (is_root_prc) CALL flinclo(fid) |
---|
4204 | ! |
---|
4205 | nbexp = 0 |
---|
4206 | ! |
---|
4207 | ! |
---|
4208 | ! Mask of permitted variables. |
---|
4209 | ! |
---|
4210 | mask(:,:) = zero |
---|
4211 | DO ip=1,iml |
---|
4212 | DO jp=1,jml |
---|
4213 | IF (soiltext(ip,jp) .GT. min_sechiba) THEN |
---|
4214 | mask(ip,jp) = un |
---|
4215 | ENDIF |
---|
4216 | ENDDO |
---|
4217 | ENDDO |
---|
4218 | ! |
---|
4219 | nbvmax = 200 |
---|
4220 | ! |
---|
4221 | callsign = "Soil types" |
---|
4222 | ! |
---|
4223 | ok_interpol = .FALSE. |
---|
4224 | DO WHILE ( .NOT. ok_interpol ) |
---|
4225 | WRITE(numout,*) "Projection arrays for ",callsign," : " |
---|
4226 | WRITE(numout,*) "nbvmax = ",nbvmax |
---|
4227 | ! |
---|
4228 | ALLOC_ERR=-1 |
---|
4229 | ALLOCATE(solt(nbvmax), STAT=ALLOC_ERR) |
---|
4230 | IF (ALLOC_ERR/=0) THEN |
---|
4231 | WRITE(numout,*) "ERROR IN ALLOCATION of solt : ",ALLOC_ERR |
---|
4232 | STOP |
---|
4233 | ENDIF |
---|
4234 | ALLOC_ERR=-1 |
---|
4235 | ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR) |
---|
4236 | IF (ALLOC_ERR/=0) THEN |
---|
4237 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR |
---|
4238 | STOP |
---|
4239 | ENDIF |
---|
4240 | sub_index(:,:,:)=0 |
---|
4241 | ALLOC_ERR=-1 |
---|
4242 | ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR) |
---|
4243 | IF (ALLOC_ERR/=0) THEN |
---|
4244 | WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR |
---|
4245 | STOP |
---|
4246 | ENDIF |
---|
4247 | sub_area(:,:)=zero |
---|
4248 | ! |
---|
4249 | CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
4250 | & iml, jml, lon_rel, lat_rel, mask, callsign, & |
---|
4251 | & nbvmax, sub_index, sub_area, ok_interpol) |
---|
4252 | ! |
---|
4253 | IF ( .NOT. ok_interpol ) THEN |
---|
4254 | DEALLOCATE(sub_area) |
---|
4255 | DEALLOCATE(sub_index) |
---|
4256 | DEALLOCATE(solt) |
---|
4257 | ! |
---|
4258 | nbvmax = nbvmax * 2 |
---|
4259 | ENDIF |
---|
4260 | ENDDO |
---|
4261 | ! |
---|
4262 | DO ib =1, nbpt |
---|
4263 | ! |
---|
4264 | soiltype(ib,:) = zero |
---|
4265 | clayfraction(ib) = zero |
---|
4266 | ! |
---|
4267 | ! GO through the point we have found |
---|
4268 | ! |
---|
4269 | ! |
---|
4270 | fopt = COUNT(sub_area(ib,:) > zero) |
---|
4271 | ! |
---|
4272 | ! Check that we found some points |
---|
4273 | ! |
---|
4274 | IF ( fopt .EQ. 0) THEN |
---|
4275 | nbexp = nbexp + 1 |
---|
4276 | soiltype(ib,:) = soiltype_default(:) |
---|
4277 | clayfraction(ib) = clayfraction_default |
---|
4278 | ELSE |
---|
4279 | ! |
---|
4280 | DO ilf = 1,fopt |
---|
4281 | solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2)) |
---|
4282 | ENDDO |
---|
4283 | ! |
---|
4284 | sgn = zero |
---|
4285 | ! |
---|
4286 | ! Compute the average bare soil albedo parameters |
---|
4287 | ! |
---|
4288 | DO ilf = 1,fopt |
---|
4289 | ! |
---|
4290 | ! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean |
---|
4291 | ! |
---|
4292 | IF ( (solt(ilf) .LE. classnb) .AND. (solt(ilf) .GT. 0) .AND.& |
---|
4293 | & (solt(ilf) .NE. 6)) THEN |
---|
4294 | SELECTCASE(solt(ilf)) |
---|
4295 | CASE(1) |
---|
4296 | soiltype(ib,1) = soiltype(ib,1) + sub_area(ib,ilf) |
---|
4297 | CASE(2) |
---|
4298 | soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) |
---|
4299 | CASE(3) |
---|
4300 | soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) |
---|
4301 | CASE(4) |
---|
4302 | soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) |
---|
4303 | CASE(5) |
---|
4304 | soiltype(ib,3) = soiltype(ib,3) + sub_area(ib,ilf) |
---|
4305 | CASE(7) |
---|
4306 | soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) |
---|
4307 | CASE DEFAULT |
---|
4308 | WRITE(numout,*) 'We should not be here, an impossible case appeared' |
---|
4309 | STOP 'slowproc_soilt' |
---|
4310 | END SELECT |
---|
4311 | clayfraction(ib) = clayfraction(ib) + & |
---|
4312 | & textfrac_table(solt(ilf),3) * sub_area(ib,ilf) |
---|
4313 | sgn = sgn + sub_area(ib,ilf) |
---|
4314 | ELSE |
---|
4315 | IF (solt(ilf) .GT. classnb) THEN |
---|
4316 | WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program' |
---|
4317 | STOP 'slowproc_soilt' |
---|
4318 | ENDIF |
---|
4319 | ENDIF |
---|
4320 | ! |
---|
4321 | ENDDO |
---|
4322 | ! |
---|
4323 | ! Normalize the surface |
---|
4324 | ! |
---|
4325 | IF ( sgn .LT. min_sechiba) THEN |
---|
4326 | nbexp = nbexp + 1 |
---|
4327 | soiltype(ib,:) = soiltype_default(:) |
---|
4328 | clayfraction(ib) = clayfraction_default |
---|
4329 | ELSE |
---|
4330 | soiltype(ib,:) = soiltype(ib,:)/sgn |
---|
4331 | clayfraction(ib) = clayfraction(ib)/sgn |
---|
4332 | ENDIF |
---|
4333 | ! |
---|
4334 | ENDIF |
---|
4335 | ! |
---|
4336 | ENDDO |
---|
4337 | ! |
---|
4338 | IF ( nbexp .GT. 0 ) THEN |
---|
4339 | WRITE(numout,*) 'slowproc_soilt : The interpolation of the bare soil albedo had ', nbexp |
---|
4340 | WRITE(numout,*) 'slowproc_soilt : points without data. This are either coastal points or' |
---|
4341 | WRITE(numout,*) 'slowproc_soilt : ice covered land.' |
---|
4342 | WRITE(numout,*) 'slowproc_soilt : The problem was solved by using the default soil types.' |
---|
4343 | ENDIF |
---|
4344 | ! |
---|
4345 | DEALLOCATE (lat_rel) |
---|
4346 | DEALLOCATE (lon_rel) |
---|
4347 | DEALLOCATE (mask) |
---|
4348 | DEALLOCATE (sub_area) |
---|
4349 | DEALLOCATE (sub_index) |
---|
4350 | DEALLOCATE (soiltext) |
---|
4351 | DEALLOCATE (solt) |
---|
4352 | ! |
---|
4353 | ! |
---|
4354 | RETURN |
---|
4355 | ! |
---|
4356 | END SUBROUTINE slowproc_soilt |
---|
4357 | ! |
---|
4358 | |
---|
4359 | |
---|
4360 | END MODULE slowproc |
---|