Context Navigation

Back to DevelopmentActivities/MergeHydro

DevelopmentActivities/MergeHydro: diff_1_9_4_1_Matthieu_version_Jan.diff

File diff_1_9_4_1_Matthieu_version_Jan.diff, 184.7 KB (added by nvuilsce, 13 years ago)
Fichier des diffs entre la version LMD et la version 1.9.4.1

Line
1	Seulement dans /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE: Doxyfile_ORCHIDEE
2	_______________________________________________________________________________________________________________________
3	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_global/AA_make /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_global/AA_make
4	2c2
5	< #- $Id: AA_make,v 1.4 2008/01/08 11:49:07 ssipsl Exp $
6	---
7	> #- $Id: AA_make,v 1.4.6.1 2009/11/10 14:16:45 ssipsl Exp $
8	40a41
9	> #-Q- sx9mercure mv $*.mod $(MODDIR)
10	_______________________________________________________________________________________________________________________
11	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_global/grid.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_global/grid.f90
12	26,27d25
13	< ! earth radius
14	< REAL(r_std), PARAMETER :: R_Earth = 6378000.
15	_______________________________________________________________________________________________________________________
16	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_global/interpol_help.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_global/interpol_help.f90
17	34,35d33
18	< REAL(r_std), PARAMETER :: R_Earth = 6378000.
19	< REAL(r_std), PARAMETER :: mincos = 0.0001
20	339a338,342
21	> IF ( jp == 0 ) THEN
22	> WRITE(,) "We have a big problem in the interpolation : ", ip, jp
23	> WRITE(,) "lat_rel, lalow_rel, laup_rel : ", lat_rel(ip,jp), lalow_rel(ip,jp), laup_rel(ip,jp)
24	> WRITE(,) "lonrel, lolowrel, louprel : ", lonrel, lolowrel, louprel
25	> ENDIF
26	357c360
27	< WRITE(numout,*) 'Reached value ', fopt,' for fopt on point', ib
28	---
29	> WRITE(numout,*) 'Reached value ', fopt,' for fopt on point', lalo(ib,2), lalo(ib,1)
30	383a387,398
31	> !
32	> DO ib=1,nbpt
33	> DO fopt=1,incmax
34	> IF (( indinc(ib,fopt,1) == 0 .AND. indinc(ib,fopt,2) > 0) .OR.&
35	> & ( indinc(ib,fopt,2) == 0 .AND. indinc(ib,fopt,1) > 0) ) THEN
36	> WRITE(,) "aggregate_2d PROBLEM : point =",ib, fopt," Indicies = ", &
37	> & indinc(ib,fopt,1), indinc(ib,fopt,2), areaoverlap(ib,fopt)
38	> ENDIF
39	> ENDDO
40	> ENDDO
41	>
42	>
43	440c455
44	< INTEGER(i_std) :: ip, ib, i, j, itmp, iprog, nbind
45	---
46	> INTEGER(i_std) :: ip, ib, i, j, itmp, iprog, nbind, pp, ipp
47	442c457,460
48	< INTEGER(i_std) :: ff(1)
49	---
50	> REAL(r_std) :: minlon, minlat, mini
51	> INTEGER(i_std) :: ff(1), incp
52	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: fine_ind
53	> INTEGER(i_std) :: pos_pnt(5)
54	444c462,463
55	<
56	---
57	> INTEGER :: ALLOC_ERR
58	> !
59	472a492,551
60	> ! Find appropriate resolution for index table
61	> !
62	> ff=MINLOC(resolution(:,1))
63	> coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )pi/180. R_Earth
64	> minlon=resolution(ff(1),1)/(2.0*coslat)
65	> ff=MINLOC(resolution(:,2))
66	> coslat = pi/180. * R_Earth
67	> minlat=resolution(ff(1),2)/(2.0*coslat)
68	> mini=MIN(minlon, minlat)
69	> !
70	> ! This interpolation only works if the model grid is coarser than the data grid
71	> !
72	> IF (MINVAL(resolution(:,1)) < resol_lon .OR. MINVAL(resolution(:,2)) < resol_lat) THEN
73	> WRITE(numout,*) " === WARNING == "
74	> WRITE(numout,*) "Resolution minima of the model (lon, lat) : ", &
75	> & MINVAL(resolution(:,1)), MINVAL(resolution(:,2))
76	> WRITE(numout,*) "Resolution of the file to be interpolated (fine grid) : ", resol_lon, resol_lat
77	> WRITE(numout,*) "This interpolation assumes that we aggregate from a fine grid to a coarser grid"
78	> WRITE(numout,*) "In the data submitted it apears that the model is runing on a finer grid than the data"
79	> ENDIF
80	> !
81	> incp = 10
82	> IF (mini < 0.1) THEN
83	> incp=100
84	> ELSE IF (mini < 0.01) THEN
85	> incp = 1000
86	> ENDIF
87	> !
88	> ! Allocate the needed memory for fine_ind
89	> !
90	> ALLOC_ERR=-1
91	> ALLOCATE (fine_ind(NINT(domain_minlonincp)-2:NINT(domain_maxlonincp)+2, &
92	> & NINT(domain_minlatincp)-2:NINT(domain_maxlatincp)+2), STAT=ALLOC_ERR)
93	> IF (ALLOC_ERR/=0) THEN
94	> WRITE(numout,*) "aggregate_vec : ERROR IN ALLOCATION of fine_ind : ",ALLOC_ERR
95	> STOP
96	> ENDIF
97	> !
98	> ! Generate a quick access table for the coarse grid
99	> !
100	> fine_ind(:,:) = zero
101	> !
102	> DO ib=1,nbpt
103	> coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )pi/180. R_Earth
104	> !
105	> lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
106	> lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
107	> !
108	> coslat = pi/180. * R_Earth
109	> !
110	> lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
111	> lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
112	> !
113	> fine_ind(NINT(lon_lowincp):NINT(lon_upincp),NINT(lat_lowincp):NINT(lat_upincp))=ib
114	> !
115	> ENDDO
116	> !
117	> WRITE(numout,*) 'Domaine LON range : ', domain_minlon, domain_maxlon
118	> WRITE(numout,*) 'Domaine LAT range : ', domain_minlat, domain_maxlat
119	> !
120	491,492c570,571
121	< louprel = lon_rel(ip) + resol_lon/(2.0*coslat)
122	< lolowrel = lon_rel(ip) - resol_lon/(2.0*coslat)
123	---
124	> louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), 180.)
125	> lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), -180.)
126	494,495c573,574
127	< lauprel = lat_rel(ip) + resol_lat/(2.0*coslat)
128	< lalowrel = lat_rel(ip) - resol_lat/(2.0*coslat)
129	---
130	> lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), 90.)
131	> lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), -90.)
132	505c584
133	< WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid'
134	---
135	> WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid and ', nbpt, 'for the coarse grid'
136	508,511d586
137	< #ifdef INTERPOL_ADVANCE
138	< WRITE(numout,'(2a40)')'0%--------------------------------------', &
139	< & '------------------------------------100%'
140	< #endif
141	518,545d592
142	< ! We select here the case which is fastest, i.e. when the smaller loop is inside
143	< !
144	< IF ( nbpt > nbind ) THEN
145	< !
146	< loopnbpt : DO ib =1, nbpt
147	< !
148	< ! Give a progress meter
149	< !
150	< #ifdef INTERPOL_ADVANCE
151	< iprog = NINT(REAL(ib,r_std)/REAL(nbpt,r_std)79.) - NINT(REAL(ib-1,r_std)/REAL(nbpt,r_std)79.)
152	< IF ( iprog .NE. 0 ) THEN
153	< WRITE(numout,'(a1,$)') 'x'
154	< ENDIF
155	< #endif
156	< !
157	< ! We find the 4 limits of the grid-box. As we transform the resolution of the model
158	< ! into longitudes and latitudes we do not have the problem of periodicity.
159	< ! coslat is a help variable here !
160	< !
161	< coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )pi/180. R_Earth
162	< !
163	< lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
164	< lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
165	< !
166	< coslat = pi/180. * R_Earth
167	< !
168	< lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
169	< lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
170	547,550d593
171	< ! Find the grid boxes from the data that go into the model's boxes
172	< ! We still work as if we had a regular grid ! Well it needs to be localy regular so
173	< ! so that the longitude at the latitude of the last found point is close to the one
174	< ! of the next point.
175	553a597
176	> !
177	562,563c606,607
178	< louprel = lon_rel(ip) + resol_lon/(2.0*coslat)
179	< lolowrel = lon_rel(ip) - resol_lon/(2.0*coslat)
180	---
181	> louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), domain_maxlon)
182	> lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), domain_minlon)
183	567,568c611,612
184	< lauprel = lat_rel(ip) + resol_lat/(2.0*coslat)
185	< lalowrel = lat_rel(ip) - resol_lat/(2.0*coslat)
186	---
187	> lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), domain_maxlat)
188	> lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), domain_minlat)
189	571,575c615,617
190	< IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
191	< & lolowrel < lon_low .AND. louprel > lon_low .OR. &
192	< & lolowrel < lon_up .AND. louprel > lon_up ) THEN
193	< !
194	< ! Now that we have the longitude let us find the latitude
195	---
196	> pos_pnt(:) = zero
197	> ipp = zero
198	> pp = fine_ind(NINT(lonrelincp),NINT(latrelincp))
199	577,602c619,621
200	< IF ( latrel > lat_low .AND. latrel < lat_up .OR. &
201	< & lalowrel < lat_low .AND. lauprel > lat_low .OR.&
202	< & lalowrel < lat_up .AND. lauprel > lat_up) THEN
203	< !
204	< fopt(ib) = fopt(ib) + 1
205	< fopt_max = MAX ( fopt(ib), fopt_max )
206	< IF ( fopt(ib) > incmax) THEN
207	< err_fopt=.TRUE.
208	< EXIT loopnbpt
209	< ELSE
210	< !
211	< ! Get the area of the fine grid in the model grid
212	< !
213	< coslat = MAX( COS( lat_rel(ip) * pi/180. ), mincos )
214	< ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))pi/180. R_Earth * coslat
215	< ay = (MIN(lat_up, lauprel)-MAX(lat_low,lalowrel))pi/180. R_Earth
216	< !
217	< areaoverlap(ib, fopt(ib)) = ax*ay
218	< indinc(ib, fopt(ib)) = ip
219	< !
220	< ! If this point was 100% within the grid then we can de-select it from our
221	< ! list as it can not be in another mesh of the coarse grid.
222	< !
223	< IF ( louprel < lon_up .AND. lolowrel > lon_low .AND. &
224	< & lauprel < lat_up .AND. lalowrel > lat_low ) THEN
225	< searchind(i) = 0
226	---
227	> IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
228	> pos_pnt(ipp+1) = pp
229	> ipp = ipp + 1
230	603a623
231	> pp = fine_ind(NINT(louprelincp),NINT(lauprelincp))
232	604a625,627
233	> IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
234	> pos_pnt(ipp+1) = pp
235	> ipp = ipp + 1
236	606,610c629
237	< ENDIF ! IF lat
238	< ENDIF ! IF lon
239	< ENDDO
240	< !
241	< ! De-select the marked points
242	---
243	> pp = fine_ind(NINT(louprelincp),NINT(lalowrelincp))
244	612,617c631,633
245	< itmp = nbind
246	< nbind = 0
247	< DO i=1,itmp
248	< IF ( searchind(i) > 0 ) THEN
249	< nbind = nbind + 1
250	< searchind(nbind) = searchind(i)
251	---
252	> IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
253	> pos_pnt(ipp+1) = pp
254	> ipp = ipp + 1
255	619c635
256	< ENDDO
257	---
258	> pp = fine_ind(NINT(lolowrelincp),NINT(lauprelincp))
259	621,633c637,639
260	< ENDDO loopnbpt
261	< IF (err_fopt) THEN
262	< WRITE(numout,*) 'Reached value ', fopt(ib),' for fopt on point', ib
263	< CALL ipslerr(2,'aggregate_vec (nbpt > nbind)', &
264	< 'Working on variable :'//callsign, &
265	< 'Reached incmax value for fopt.',&
266	< 'Please increase incmax in subroutine calling aggregate')
267	< IF (PRESENT(ok)) THEN
268	< ok = .FALSE.
269	< RETURN
270	< ELSE
271	< STOP "aggregate_vec"
272	< ENDIF
273	---
274	> IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
275	> pos_pnt(ipp+1) = pp
276	> ipp = ipp + 1
277	634a641
278	> pp = fine_ind(NINT(lolowrelincp),NINT(lalowrelincp))
279	636,648c643,645
280	< ELSE
281	< !
282	< ib = 1
283	< !
284	< loopnbind : DO i=1,nbind
285	< !
286	< !
287	< ! Give a progress meter
288	< !
289	< #ifdef INTERPOL_ADVANCE
290	< iprog = NINT(REAL(i,r_std)/REAL(nbind,r_std)79.) - NINT(REAL(i-1,r_std)/REAL(nbind,r_std)79.)
291	< IF ( iprog .NE. 0 ) THEN
292	< WRITE(numout,'(a1,$)') 'y'
293	---
294	> IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
295	> pos_pnt(ipp+1) = pp
296	> ipp = ipp + 1
297	650,661d646
298	< #endif
299	< !
300	< ip = searchind(i)
301	< !
302	< ! Either the center of the data grid point is in the interval of the model grid or
303	< ! the East and West limits of the data grid point are on either sides of the border of
304	< ! the data grid.
305	< !
306	< lonrel = lon_rel(ip)
307	< coslat = MAX(COS(lat_rel(ip) * pi/180. ), mincos )pi/180. R_Earth
308	< louprel = lon_rel(ip) + resol_lon/(2.0*coslat)
309	< lolowrel = lon_rel(ip) - resol_lon/(2.0*coslat)
310	663,666d647
311	< latrel = lat_rel(ip)
312	< coslat = pi/180. * R_Earth
313	< lauprel = lat_rel(ip) + resol_lat/(2.0*coslat)
314	< lalowrel = lat_rel(ip) - resol_lat/(2.0*coslat)
315	667a649
316	> IF ( ipp > zero ) THEN
317	669,674c651,652
318	< found = .FALSE.
319	< j = 1
320	< !
321	< DO WHILE ( .NOT. found .AND. j <= nbpt )
322	< ! Just count the number of time we were through
323	< j = j+1
324	---
325	> DO pp=1,ipp
326	> ib = pos_pnt(pp)
327	689a668,669
328	> !
329	> !
330	701a682
331	> !
332	704d684
333	< EXIT loopnbind
334	715d694
335	< found = .TRUE.
336	718,721c697,701
337	< ENDIF ! IF lat
338	< ENDIF ! IF lon
339	< ENDDO ! While loop
340	< ENDDO loopnbind
341	---
342	> ENDIF
343	> ENDIF
344	> ENDDO
345	> ENDIF
346	> ENDDO
347	737,740d716
348	< IF ( .NOT. found ) THEN
349	< ! We need to step on in the coarse grid
350	< ib = ib + 1
351	< IF ( ib > nbpt ) ib = ib-nbpt
352	742,743d717
353	< ENDIF
354	< ENDIF
355	745c719,723
356	< WRITE(numout,*) "aggregate_vec nbvmax = ",incmax, "max used = ",fopt_max
357	---
358	> WRITE(numout,*) "aggregate_vec : ",COUNT(fopt(:) .EQ. zero), &
359	> & "did not find any corresponding data in the input file."
360	> WRITE(numout,) "aggregate_vec : This is ", COUNT(fopt(:) .EQ. zero)/FLOAT(nbpt)100., &
361	> & " % of the grid"
362	> WRITE(numout,*) "aggregate_vec : nbvmax = ",incmax, "max used = ",fopt_max
363	749a728
364	> DEALLOCATE (fine_ind)
365	817a797,798
366	> INTEGER(i_std) :: ib, iv
367	>
368	820a802,805
369	>
370	> sub_index(:,:,:) = 0
371	> sub_area(:,:) = zero
372	>
373	825d809
374	<
375	_______________________________________________________________________________________________________________________
376	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parallel/AA_make /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parallel/AA_make
377	2c2
378	< #- $Id: AA_make,v 1.5 2008/01/08 11:49:07 ssipsl Exp $
379	---
380	> #- $Id: AA_make,v 1.5.6.1 2009/11/10 14:16:45 ssipsl Exp $
381	40a41
382	> #-Q- sx9mercure mv $*.mod $(MODDIR)
383	_______________________________________________________________________________________________________________________
384	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parallel/data_para.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parallel/data_para.f90
385	128a129,145
386	> ELSE IF (MPI_VERSION == 2) THEN
387	> ! Version MPI 2.x
388	> IF (i_std==i_4) THEN
389	> MPI_INT_ORCH=MPI_INTEGER
390	> ELSEIF (i_std==i_8) THEN
391	> MPI_INT_ORCH=MPI_INTEGER
392	> ELSE
393	> MPI_INT_ORCH=MPI_INTEGER
394	> ENDIF
395	>
396	> IF (r_std==r_4) THEN
397	> MPI_REAL_ORCH=MPI_REAL
398	> ELSEIF (r_std==r_8) THEN
399	> MPI_REAL_ORCH=MPI_DOUBLE_PRECISION
400	> ELSE
401	> MPI_REAL_ORCH=MPI_REAL
402	> ENDIF
403	_______________________________________________________________________________________________________________________
404	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parallel/transfert_para.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parallel/transfert_para.f90
405	569c569
406	< LOGICAL, PARAMETER :: check=.FALSE.
407	---
408	> LOGICAL, PARAMETER :: check=.TRUE.
409	_______________________________________________________________________________________________________________________
410	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parameters/AA_make /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parameters/AA_make
411	2c2
412	< #- $Id: AA_make,v 1.16 2008/01/08 11:49:07 ssipsl Exp $
413	---
414	> #- $Id: AA_make,v 1.16.6.1 2009/11/10 14:16:45 ssipsl Exp $
415	42a43
416	> #-Q- sx9mercure mv $*.mod $(MODDIR)
417	_______________________________________________________________________________________________________________________
418	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parameters/constantes.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parameters/constantes.f90
419	37a38,40
420	> REAL(r_std), PARAMETER :: R_Earth = 6378000.
421	> REAL(r_std), PARAMETER :: mincos = 0.0001
422	> !
423	44a48
424	> INTEGER(i_std),PARAMETER :: undef_int = 999999999
425	71a76,77
426	> LOGICAL :: do_floodplains
427	> LOGICAL :: do_irrigation
428	_______________________________________________________________________________________________________________________
429	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parameters/constantes_soil.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parameters/constantes_soil.f90
430	11a12,15
431	> !
432	> LOGICAL, SAVE :: check_waterbal=.TRUE. !! The check the water balance
433	> LOGICAL, SAVE :: waterbal_error=.FALSE. !! If true the water balance is bad
434	> !
435	23a28,39
436	> ! Maximum depth of soil reservoir. It is fixed accordingly to nslm above (11 -> 2)
437	> ! If a depth map is given, nbdl, and nslm will be put to nslm_max and dpu_max will be changed
438	> ! in intersurf
439	> REAL(r_std),SAVE :: dpu_max=2.0_r_std
440	> ! Number of levels min and max in CWRR (used only when a depth map is given)
441	> INTEGER(i_std),PARAMETER :: nslm_min=10
442	> INTEGER(i_std),PARAMETER :: nslm_max=13
443	> !
444	> !Number of soil classes
445	> !
446	> ! Number of soil textures (Silt, Sand, Clay)
447	> INTEGER(i_std),PARAMETER :: ntext=3
448	25a42,47
449	> !
450	> ! For FAO Classification
451	> INTEGER(i_std),PARAMETER :: nscm_fao=3
452	> ! For USDA Classification
453	> INTEGER(i_std),PARAMETER :: nscm_usda=12
454	> INTEGER(i_std), SAVE :: nscm=nscm_fao
455	90a113,126
456	> !-
457	> !- 1. Parameters for FAO Classification
458	> !-
459	>
460	> !-
461	> ! Parameters for soil type distribution
462	> !-
463	> ! Default soil texture distribution in the following order :
464	> ! COARSE, MEDIUM, FINE
465	> REAL(r_std),DIMENSION(nscm_fao),SAVE :: soilclass_default_fao = &
466	> & (/ 0.28, 0.52, 0.20 /)
467	> REAL(r_std),DIMENSION(nscm_fao*2-1),SAVE :: soilclass_default_fao2 = &
468	> & (/ 0.28, 0.0, 0.52, 0.0, 0.20 /)
469	> INTEGER, SAVE :: jsc_default = 2
470	92c128
471	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: nvan = &
472	---
473	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: nvan_fao = &
474	94,99c130,131
475	< !!$! Van genuchten coefficient a (cm^{-1})
476	< !!$ REAL(r_std),PARAMETER,DIMENSION(nstm) :: avan = &
477	< !!$ & (/ 0.036_r_std, 0.036_r_std, 0.036_r_std /)
478	< !TdO
479	< ! Van genuchten coefficient a (mm^{-1})
480	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: avan = &
481	---
482	> ! Van genuchten coefficient a (cm^{-1}) BIG BUG -> mm^{-1}
483	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: avan_fao = &
484	101,106c133,134
485	< ! CWRR linearisation
486	< INTEGER(i_std),PARAMETER :: imin = 1
487	< ! number of interval for CWRR
488	< INTEGER(i_std),PARAMETER :: nbint = 100
489	< ! number of points for CWRR
490	< INTEGER(i_std),PARAMETER :: imax = nbint+1
491	---
492	> ! & (/ 0.075_r_std, 0.036_r_std, 0.019_r_std /)
493	> ! & (/ 0.036_r_std, 0.036_r_std, 0.036_r_std /)
494	108c136
495	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcr = &
496	---
497	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcr_fao = &
498	111c139
499	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcs = &
500	---
501	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcs_fao = &
502	113,116d140
503	< ! Total depth of soil reservoir (m)
504	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: dpu = &
505	< & (/ 2.0_r_std, 2.0_r_std, 2.0_r_std /)
506	< ! dpu must be constant over the different soil types
507	118c142
508	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: ks = &
509	---
510	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: ks_fao = &
511	120,121c144,145
512	< ! Soil moisture above which transpir is max
513	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: pcent = &
514	---
515	> ! Fraction of saturated volumetric soil moisture above which transpir is max
516	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: pcent_fao = &
517	124c148
518	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: free_drain_max = &
519	---
520	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: free_drain_max_fao = &
521	127c151
522	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcf = &
523	---
524	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcf_fao = &
525	130c154
526	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcw = &
527	---
528	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcw_fao = &
529	131a156
530	> ! & (/ 0.07_r_std, 0.085_r_std, 0.10_r_std /)
531	133c158
532	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mc_awet = &
533	---
534	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mc_awet_fao = &
535	136c161
536	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: mc_adry = &
537	---
538	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mc_adry_fao = &
539	139c164
540	< REAL(r_std),PARAMETER,DIMENSION(nstm) :: psis = &
541	---
542	> REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: psis_fao = &
543	140a166,248
544	> !-
545	> !- 2. Parameters for USDA Classification
546	> !-
547	>
548	> !-
549	> ! Parameters for soil type distribution
550	> !-
551	> ! Default soil texture distribution in the following order :
552	> ! sand, loam and clay ??? OR COARSE, MEDIUM, FINE???
553	> REAL(r_std),DIMENSION(nscm_usda),SAVE :: soilclass_default_usda = &
554	> & (/ 0.28, 0.52, 0.20, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
555	>
556	> ! Van genuchten coefficient n
557	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: nvan_usda = &
558	> & (/ 2.68_r_std, 2.28_r_std, 1.89_r_std, 1.41_r_std, &
559	> & 1.37_r_std, 1.56_r_std, 1.48_r_std, 1.23_r_std, &
560	> & 1.31_r_std, 1.23_r_std, 1.09_r_std, 1.09_r_std /)
561	> ! Van genuchten coefficient a (cm^{-1}) BIG BUG!!! -> mm^{-1}
562	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: avan_usda = &
563	> & (/ 0.0145_r_std, 0.0124_r_std, 0.0075_r_std, 0.0020_r_std, &
564	> & 0.0016_r_std, 0.0036_r_std, 0.0059_r_std, 0.0010_r_std, &
565	> & 0.0019_r_std, 0.0027_r_std, 0.0005_r_std, 0.0008_r_std /)
566	> ! & (/ 0.145_r_std, 0.124_r_std, 0.075_r_std, 0.020_r_std, &
567	> ! & 0.016_r_std, 0.036_r_std, 0.059_r_std, 0.010_r_std, &
568	> ! & 0.019_r_std, 0.027_r_std, 0.005_r_std, 0.008_r_std /)
569	> ! Sand, Loamy Sand, Sandy Loam, Silt Loam, Silt, Loam, Sandy Clay Loam, Silty Clay Loam, Clay Loam, Sandy Clay, Silty Clay, Clay
570	> ! Residual soil water content
571	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcr_usda = &
572	> & (/ 0.045_r_std, 0.057_r_std, 0.065_r_std, 0.067_r_std, &
573	> & 0.034_r_std, 0.078_r_std, 0.100_r_std, 0.089_r_std, &
574	> & 0.095_r_std, 0.100_r_std, 0.070_r_std, 0.068_r_std /)
575	> ! Saturated soil water content
576	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcs_usda = &
577	> & (/ 0.43_r_std, 0.41_r_std, 0.41_r_std, 0.45_r_std, &
578	> & 0.46_r_std, 0.43_r_std, 0.39_r_std, 0.43_r_std, &
579	> & 0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std /)
580	> ! Hydraulic conductivity Saturation (mm/d)
581	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: ks_usda = &
582	> & (/ 7128.0_r_std, 3501.6_r_std, 1060.8_r_std, 108.0_r_std, &
583	> & 60.0_r_std, 249.6_r_std, 314.4_r_std, 16.8_r_std, &
584	> & 62.4_r_std, 28.8_r_std, 4.8_r_std, 48.0_r_std /)
585	> ! Soil moisture above which transpir is max
586	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: pcent_usda = &
587	> & (/ 0.5_r_std, 0.5_r_std, 0.5_r_std, 0.5_r_std, &
588	> & 0.5_r_std, 0.5_r_std, 0.5_r_std, 0.5_r_std, &
589	> & 0.5_r_std, 0.5_r_std, 0.5_r_std, 0.5_r_std /)
590	> ! Max value of the permeability coeff at the bottom of the soil
591	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: free_drain_max_usda = &
592	> & (/ 1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std, &
593	> & 1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std, &
594	> & 1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std /)
595	> ! Volumetric water content field capacity
596	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcf_usda = &
597	> & (/ 0.32_r_std, 0.32_r_std, 0.32_r_std, 0.32_r_std, &
598	> & 0.32_r_std, 0.32_r_std, 0.32_r_std, 0.32_r_std, &
599	> & 0.32_r_std, 0.32_r_std, 0.32_r_std, 0.32_r_std /)
600	> ! Volumetric water content Wilting pt
601	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcw_usda = &
602	> & (/ 0.10_r_std, 0.10_r_std, 0.10_r_std, 0.10_r_std, &
603	> & 0.10_r_std, 0.10_r_std, 0.10_r_std, 0.10_r_std, &
604	> & 0.10_r_std, 0.10_r_std, 0.10_r_std, 0.10_r_std /)
605	> ! Vol. wat. cont. above which albedo is cst
606	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mc_awet_usda = &
607	> & (/ 0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std, &
608	> & 0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std, &
609	> & 0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std /)
610	> ! Vol. wat. cont. below which albedo is cst
611	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mc_adry_usda = &
612	> & (/ 0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &
613	> & 0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &
614	> & 0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std /)
615	> ! Matrix potential at saturation (mm)
616	> REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: psis_usda = &
617	> & (/ -300.0_r_std, -300.0_r_std, -300.0_r_std, -300.0_r_std, &
618	> & -300.0_r_std, -300.0_r_std, -300.0_r_std, -300.0_r_std, &
619	> & -300.0_r_std, -300.0_r_std, -300.0_r_std, -300.0_r_std /)
620	>
621	> ! CWRR linearisation
622	> INTEGER(i_std),PARAMETER :: imin = 1
623	> ! number of interval for CWRR
624	> INTEGER(i_std),PARAMETER :: nbint = 50
625	> ! number of points for CWRR
626	> INTEGER(i_std),PARAMETER :: imax = nbint+1
627	_______________________________________________________________________________________________________________________
628	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_parameters/constantes_veg.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_parameters/constantes_veg.f90
629	173,177d172
630	< !!$ CHARACTER(len=5),DIMENSION(nvm),SAVE :: type_of_lai = &
631	< !!$ & (/ 'mean ', 'mean ', 'inter', 'mean ', 'mean ', &
632	< !!$ & 'inter', 'mean ', 'inter', 'inter', 'inter', &
633	< !!$ & 'inter', 'inter', 'inter' /)
634	< ! Test Nathalie : Even Sempervirens vegetation is allowed to have a small seasonal cycle.
635	179,180c174,175
636	< & (/ 'inter', 'inter', 'inter', 'inter', 'inter', &
637	< & 'inter', 'inter', 'inter', 'inter', 'inter', &
638	---
639	> & (/ 'mean ', 'mean ', 'inter', 'mean ', 'mean ', &
640	> & 'inter', 'mean ', 'inter', 'inter', 'inter', &
641	233c228
642	< INTEGER(i_std),DIMENSION(nvm,nstm) :: pref_soil_veg
643	---
644	> INTEGER(i_std),DIMENSION(nvm) :: pref_soil_veg
645	_______________________________________________________________________________________________________________________
646	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_sechiba/AA_make /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba/AA_make
647	2c2
648	< #- $Id: AA_make,v 1.20 2008/01/08 11:49:07 ssipsl Exp $
649	---
650	> #- $Id: AA_make,v 1.20.6.1 2009/11/10 14:16:45 ssipsl Exp $
651	84a85
652	> #-Q- sx9mercure mv $*.mod $(MODDIR)
653	_______________________________________________________________________________________________________________________
654	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_sechiba/intersurf.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba/intersurf.f90
655	1d0
656	<
657	38c37
658	< MODULE PROCEDURE intersurf_main_2d, intersurf_main_1d, intersurf_gathered, intersurf_gathered_2m
659	---
660	> MODULE PROCEDURE intersurf_main_2d, intersurf_main_1d, intersurf_gathered
661	185a185,187
662	> LOGICAL, SAVE :: fatmco2 !! Flag to force the value of atmospheric CO2 for vegetation.
663	> REAL(r_std), SAVE :: atmco2 !! atmospheric CO2
664	> !
665	225a228,254
666	> !Config Key = FORCE_CO2_VEG
667	> !Config Desc = Flag to force the value of atmospheric CO2 for vegetation.
668	> !Config Def = FALSE
669	> !Config Help = If this flag is set to true, the ATM_CO2 parameter is used
670	> !Config to prescribe the atmospheric CO2.
671	> !Config This Flag is only use in couple mode.
672	> !
673	> fatmco2=.FALSE.
674	> CALL getin_p('FORCE_CO2_VEG',fatmco2)
675	> !
676	> ! Next flag is only use in couple mode with a gcm in intersurf.
677	> ! In forced mode, it has already been read and set in driver.
678	> IF ( fatmco2 ) THEN
679	> !Config Key = ATM_CO2
680	> !Config IF = FORCE_CO2_VEG (in not forced mode)
681	> !Config Desc = Value for atm CO2
682	> !Config Def = 350.
683	> !Config Help = Value to prescribe the atm CO2.
684	> !Config For pre-industrial simulations, the value is 286.2 .
685	> !Config 348. for 1990 year.
686	> !
687	> atmco2=350.
688	> CALL getin_p('ATM_CO2',atmco2)
689	> WRITE(numout,*) 'atmco2 ',atmco2
690	> ENDIF
691	> !
692	> !
693	287d315
694	< zccanopy(ik) = ccanopy(i,j)
695	296a325,335
696	> IF ( fatmco2 ) THEN
697	> zccanopy(:) = atmco2
698	> WRITE (numout,*) 'Modification of the ccanopy value. CO2 = ',atmco2
699	> ELSE
700	> DO ik=1, kjpindex
701	> j = ((kindex(ik)-1)/iim) + 1
702	> i = (kindex(ik) - (j-1)*iim)
703	> zccanopy(:) = ccanopy(i,j)
704	> ENDDO
705	> ENDIF
706	> !
707	362c401
708	< IF ( check ) WRITE(numout,*) 'Calling sechiba'
709	---
710	> IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_main_2d'
711	364a404
712	> !
713	370c410
714	< CALL histsync(hist_id)
715	---
716	> !!$ CALL histsync(hist_id)
717	375c415
718	< CALL histsync(hist2_id)
719	---
720	> !!$ CALL histsync(hist2_id)
721	378a419
722	> !
723	383,385c424
724	< ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul rveget
725	< ! & zzlev, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
726	< & zzlev, zu, zv, zqair, zqair, ztemp_air, ztemp_air, zepot_air, zccanopy, &
727	---
728	> & zzlev, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
729	467,469c506
730	< ! Ajout Nathalie - Juin 2006 - on conserve q2m/t2m
731	< CALL histwrite (hist_id, 'q2m', itau_sechiba, qair, iim*jjm, kindex)
732	< CALL histwrite (hist_id, 't2m', itau_sechiba, temp_air, iim*jjm, kindex)
733	---
734	> !
735	486,487d522
736	< CALL histwrite (hist2_id, 'q2m', itau_sechiba, qair, iim*jjm, kindex)
737	< CALL histwrite (hist2_id, 't2m', itau_sechiba, temp_air, iim*jjm, kindex)
738	495a531,532
739	> CALL histwrite (hist_id, 'Tair', itau_sechiba, temp_air, iim*jjm, kindex)
740	> CALL histwrite (hist_id, 'Qair', itau_sechiba, qair, iim*jjm, kindex)
741	506,508c543,545
742	< IF (dw .EQ. xrdt) THEN
743	< CALL histsync()
744	< ENDIF
745	---
746	> !!$ IF (dw .EQ. xrdt) THEN
747	> !!$ CALL histsync()
748	> !!$ ENDIF
749	640a678,680
750	> LOGICAL, SAVE :: fatmco2 !! Flag to force the value of atmospheric CO2 for vegetation.
751	> REAL(r_std), SAVE :: atmco2 !! atmospheric CO2
752	> !
753	695a736,761
754	> !Config Key = FORCE_CO2_VEG
755	> !Config Desc = Flag to force the value of atmospheric CO2 for vegetation.
756	> !Config Def = FALSE
757	> !Config Help = If this flag is set to true, the ATM_CO2 parameter is used
758	> !Config to prescribe the atmospheric CO2.
759	> !Config This Flag is only use in couple mode.
760	> !
761	> fatmco2=.FALSE.
762	> CALL getin_p('FORCE_CO2_VEG',fatmco2)
763	> !
764	> ! Next flag is only use in couple mode with a gcm in intersurf.
765	> ! In forced mode, it has already been read and set in driver.
766	> IF ( fatmco2 ) THEN
767	> !Config Key = ATM_CO2
768	> !Config IF = FORCE_CO2_VEG (in not forced mode)
769	> !Config Desc = Value for atm CO2
770	> !Config Def = 350.
771	> !Config Help = Value to prescribe the atm CO2.
772	> !Config For pre-industrial simulations, the value is 286.2 .
773	> !Config 348. for 1990 year.
774	> !
775	> atmco2=350.
776	> CALL getin_p('ATM_CO2',atmco2)
777	> WRITE(numout,*) 'atmco2 ',atmco2
778	> ENDIF
779	> !
780	733d798
781	< zccanopy(ik) = ccanopy(kindex(ik))
782	742a808,815
783	> IF ( fatmco2 ) THEN
784	> zccanopy(:) = atmco2
785	> WRITE (numout,*) 'Modification of the ccanopy value. CO2 = ',atmco2
786	> ELSE
787	> DO ik=1, kjpindex
788	> zccanopy(ik) = ccanopy(kindex(ik))
789	> ENDDO
790	> ENDIF
791	747c820
792	< IF ( check ) WRITE(numout,*) 'Calling sechiba'
793	---
794	> IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_main_1d'
795	755c828
796	< CALL histwrite(hist_id, 'LandPoints', itau_sechiba+1, (/ REAL(kindex) /), kjpindex, kindex)
797	---
798	> CALL histwrite(hist_id, 'LandPoints', itau_sechiba+1, REAL(kindex), kjpindex, kindex)
799	760c833
800	< CALL histsync(hist_id)
801	---
802	> !!$ CALL histsync(hist_id)
803	763c836
804	< CALL histwrite(hist2_id, 'LandPoints', itau_sechiba+1, (/ REAL(kindex) /), kjpindex, kindex)
805	---
806	> CALL histwrite(hist2_id, 'LandPoints', itau_sechiba+1, REAL(kindex), kjpindex, kindex)
807	765c838
808	< CALL histsync(hist2_id)
809	---
810	> !!$ CALL histsync(hist2_id)
811	773,775c846
812	< ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul rveget
813	< ! & zzlev, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
814	< & zzlev, zu, zv, zqair, zqair, ztemp_air, ztemp_air, zepot_air, zccanopy, &
815	---
816	> & zzlev, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
817	856,858c927
818	< ! Ajouts Nathalie - Juin 2006 - sauvegarde de t2m et q2m
819	< CALL histwrite (hist_id, 'q2m', itau_sechiba, qair, iim*jjm, kindex)
820	< CALL histwrite (hist_id, 't2m', itau_sechiba, temp_air, iim*jjm, kindex)
821	---
822	> !
823	875,876d943
824	< CALL histwrite (hist2_id, 'q2m', itau_sechiba, qair, iim*jjm, kindex)
825	< CALL histwrite (hist2_id, 't2m', itau_sechiba, temp_air, iim*jjm, kindex)
826	884a952,953
827	> CALL histwrite (hist_id, 'Tair', itau_sechiba, temp_air, iim*jjm, kindex)
828	> CALL histwrite (hist_id, 'Qair', itau_sechiba, qair, iim*jjm, kindex)
829	1233d1301
830	<
831	1323c1391
832	< IF ( check ) WRITE(numout,*) 'Calling sechiba'
833	---
834	> IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_gathered'
835	1370c1438
836	< CALL histsync(hist_id)
837	---
838	> !!$ CALL histsync(hist_id)
839	1375c1443
840	< CALL histsync(hist2_id)
841	---
842	> !!$ CALL histsync(hist2_id)
843	1383,1385c1451
844	< ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul rveget
845	< ! & zlev, u, v, qair, temp_air, epot_air, ccanopy, &
846	< & zlev, u, v, qair, qair, temp_air, temp_air, epot_air, zccanopy, &
847	---
848	> & zlev, u, v, qair, temp_air, epot_air, zccanopy, &
849	1486,1487d1551
850	< CALL histwrite (hist_id, 't2m', itau_sechiba, dtair, iim*jjm, kindex)
851	< CALL histwrite (hist_id, 'q2m', itau_sechiba, dqair, iim*jjm, kindex)
852	1505,1506d1568
853	< CALL histwrite (hist2_id, 't2m', itau_sechiba, dtair, iim*jjm, kindex)
854	< CALL histwrite (hist2_id, 'q2m', itau_sechiba, dqair, iim*jjm, kindex)
855	1515a1578,1579
856	> CALL histwrite (hist_id, 'Tair', itau_sechiba, dtair, iim*jjm, kindex)
857	> CALL histwrite (hist_id, 'Qair', itau_sechiba, dqair, iim*jjm, kindex)
858	1675c1739,1741
859	<
860	---
861	> ! LOCAL
862	> CHARACTER(LEN=30) :: classif
863	> !
864	1685a1752,1775
865	> !Config Key = SOILTYPE_CLASSIF
866	> !Config Desc = Type of classification used for the map of soil types
867	> !Config Def = zobler
868	> !Config If = !IMPOSE_VEG
869	> !Config Help = The classification used in the file that we use here
870	> !Config = There are three classification supported:
871	> !Config = FAO (3 soil types), Zobler (7 converted to 3) and USDA (12)
872	> !
873	> !-tdo- Suivant le type de classification utilisee pour le sol, on adapte nscm
874	> classif = 'zobler'
875	> CALL getin('SOILTYPE_CLASSIF',classif)
876	> SELECTCASE (classif)
877	> CASE ('zobler', 'fao','none')
878	> nscm = nscm_fao
879	> CASE ('fao2')
880	> nscm = 2 * nscm_fao-1
881	> CASE ('usda')
882	> nscm = nscm_usda
883	> CASE DEFAULT
884	> PRINT *, "Unsupported soil type classification. Choose between zobler, fao and usda according to the map"
885	> STOP 'intsurf_config'
886	> ENDSELECT
887	> !
888	> !
889	1743,1748c1833,1868
890	< IF ( control_flags%hydrol_cwrr ) then
891	< CALL ipslerr (2,'intsurf_config', &
892	< & 'You will use in this run the second version of CWRR hydrology in ORCHIDEE.',&
893	< & 'This model hasn''t been tested for global run yet.', &
894	< & '(check your parameters)')
895	< ENDIF
896	---
897	> WRITE(numout,*) "CWRR hydrology is activated : ",control_flags%hydrol_cwrr
898	> !
899	> !Config Key = DO_IRRIGATION
900	> !Config Desc = Should we compute an irrigation flux
901	> !Config Def = FALSE
902	> !Config Help = This parameters allows the user to ask the model
903	> !Config to compute an irigation flux. This performed for the
904	> !Config on very simple hypothesis. The idea is to have a good
905	> !Config map of irrigated areas and a simple function which estimates
906	> !Config the need to irrigate.
907	> !
908	> control_flags%do_irrigation = .FALSE.
909	> CALL getin_p('DO_IRRIGATION', control_flags%do_irrigation)
910	> !
911	> !Config Key = DO_FLOODPLAINS
912	> !Config Desc = Should we include floodplains
913	> !Config Def = FALSE
914	> !Config Help = This parameters allows the user to ask the model
915	> !Config to take into account the flood plains and return
916	> !Config the water into the soil moisture. It then can go
917	> !Config back to the atmopshere. This tried to simulate
918	> !Config internal deltas of rivers.
919	> !
920	> control_flags%do_floodplains = .FALSE.
921	> CALL getin_p('DO_FLOODPLAINS', control_flags%do_floodplains)
922	> !
923	> !Config Key = CHECK_WATERBAL
924	> !Config Desc = Should we check the global water balance
925	> !Config Def = TRUE
926	> !Config Help = This parameters allows the user to check
927	> !Config the integrated water balance at the end
928	> !Config of each time step
929	> !
930	> check_waterbal = .FALSE.
931	> CALL getin_p('CHECK_WATERBAL', check_waterbal)
932	> !
933	1890a2011
934	> !
935	2031c2152
936	< INTEGER(i_std) :: ier
937	---
938	> INTEGER(i_std) :: ier, jv
939	2104a2226,2230
940	> DO jv = 1, nslm-1
941	> diaglev(jv) = dpu_max/(2*(nslm-1) -1) ( ( 2(jv-1) -1) + ( 2(jv) -1) ) / deux
942	> ENDDO
943	> diaglev(nslm) = dpu_max
944	> !
945	2241c2367
946	< & nslm, solay, solayax_id)
947	---
948	> & nslm, diaglev(1:nslm), solayax_id)
949	2261a2388,2399
950	> !
951	> IF ( control_flags%hydrol_cwrr ) THEN
952	> CALL histdef(hist_id, 'reinf_slope', 'Slope index for each grid box', '-', &
953	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(5), dt,dw)
954	> CALL histdef(hist_id, 'soiltile', 'Fraction of soil tiles', '%', &
955	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, once(5), dt,dw)
956	> CALL histdef(hist_id, 'soilindex', 'Soil index', '-', &
957	> & iim,jjm, hori_id, 1, 1, 1, -99, 32, once(5), dt,dw)
958	> CALL histdef(hist_id, 'soildepth', 'Soil Depth', 'm', &
959	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, once(1), dt,dw)
960	> ENDIF
961	> !
962	2289,2297c2427,2431
963	< IF ( control_flags%hydrol_cwrr ) THEN
964	< CALL histdef(hist_id, 'evapnu_soil', 'Bare soil evap for soil type', 'mm/d', &
965	< & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(2), dt,dw)
966	< CALL histdef(hist_id, 'drainage_soil', 'Drainage for soil type', 'mm/d', &
967	< & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(2), dt,dw)
968	< CALL histdef(hist_id, 'transpir_soil', 'Transpir for soil type', 'mm/d', &
969	< & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(2), dt,dw)
970	< CALL histdef(hist_id, 'runoff_soil', 'Runoff for soil type', 'mm/d', &
971	< & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(2), dt,dw)
972	---
973	> IF ( control_flags%do_floodplains ) THEN
974	> CALL histdef(hist_id, 'floodout', 'Flow out of floodplains', 'mm/d', &
975	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(8), dt,dw)
976	> CALL histdef(hist_id, 'evapflo', 'Floodplains evaporation', 'mm/d', &
977	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(8), dt,dw)
978	2356a2491,2494
979	> IF ( control_flags%do_floodplains ) THEN
980	> CALL histdef(hist_id, 'flood_frac', 'Flooded fraction', '-', &
981	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
982	> ENDIF
983	2369a2508,2512
984	> ! var_name= "kfact_root_1" ... "kfact_root_3"
985	> WRITE (var_name,"('kfactroot_',i1)") jst
986	> CALL histdef(hist_id, var_name, 'Root fraction profile for soil type', '%', &
987	> & iim,jjm, hori_id, nslm, 1, nslm, solayax_id, 32, avescatter(1), dt,dw)
988	>
989	2371a2515,2517
990	> !
991	> CALL histdef(hist_id, 'frac_bare', 'Bare soil fraction for each tile', '-', &
992	> & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(5), dt,dw)
993	2386a2533,2536
994	> CALL histdef(hist_id, 'SWI', 'Soil wetness index','-', &
995	> & iim,jjm, hori_id, 1, 1, 1, -99, 32, avescatter(1), dt,dw)
996	> CALL histdef(hist_id, 'njsc', 'Soil class used for hydrology', '-', &
997	> & iim,jjm, hori_id, 1, 1, 1, -99, 32, once(1), dt,dw)
998	2443a2594,2601
999	> CALL histdef(hist_id, 'pondr', 'Ponds reservoir', 'kg/m^2', &
1000	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
1001	> CALL histdef(hist_id, 'returnflow', 'Returnflow to bottom layer', 'mm/d', &
1002	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw)
1003	> CALL histdef(hist_id, 'swampmap', 'Map of swamps', 'm^2', &
1004	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(8), dt,dw)
1005	> !
1006	> IF ( control_flags%do_irrigation ) THEN
1007	2445c2603
1008	< & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(7), dt,dw)
1009	---
1010	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(8), dt,dw)
1011	2447,2449c2605,2617
1012	< & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(7), dt,dw)
1013	< CALL histdef(hist_id, 'irrigmap', 'Map of irrigated areas', 'm^2', &
1014	< & iim,jjm, hori_id, 1,1,1, -99, 32, once(7), dt,dw)
1015	---
1016	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(8), dt,dw)
1017	> CALL histdef(hist_id, 'irrigmap', 'Map of irrigated surfaces', 'm^2', &
1018	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(8), dt,dw)
1019	> ENDIF
1020	> IF ( control_flags%do_floodplains ) THEN
1021	> CALL histdef(hist_id, 'floodmap', 'Map of floodplains', 'm^2', &
1022	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(8), dt,dw)
1023	> CALL histdef(hist_id, 'floodh', 'Floodplains height', 'mm', &
1024	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
1025	> CALL histdef(hist_id, 'floodr', 'Floodplains reservoir', 'kg/m^2', &
1026	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
1027	> ENDIF
1028	> !
1029	2464c2632
1030	< !MM
1031	---
1032	> !
1033	2470a2639,2640
1034	> CALL histdef(hist_id, 'vbeta5', 'Beta for bare soil', '-', &
1035	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
1036	2473,2474d2642
1037	< CALL histdef(hist_id, 'soiltype', 'Fraction of soil textures', '%', &
1038	< & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, once(8), dt,dw)
1039	2476a2645,2655
1040	> !
1041	> IF ( control_flags%hydrol_cwrr ) THEN
1042	> CALL histdef(hist_id, 'evapnu_soil', 'Bare soil evap for soil type', 'mm/d', &
1043	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(1), dt,dw)
1044	> CALL histdef(hist_id, 'drainage_soil', 'Drainage for soil type', 'mm/d', &
1045	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(1), dt,dw)
1046	> CALL histdef(hist_id, 'transpir_soil', 'Transpir for soil type', 'mm/d', &
1047	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(1), dt,dw)
1048	> CALL histdef(hist_id, 'runoff_soil', 'Runoff for soil type', 'mm/d', &
1049	> & iim,jjm, hori_id, nstm, 1, nstm, soltax_id, 32, fluxop(1), dt,dw)
1050	> ENDIF
1051	2528c2707
1052	< & nslm, solay, solayax_id)
1053	---
1054	> & nslm, diaglev(1:nslm), solayax_id)
1055	2582c2761
1056	< CALL histdef(hist_id, 'DelSWE', 'Change in SWE','kg/m^2',&
1057	---
1058	> CALL histdef(hist_id, 'DelSurfStor', 'Change in Surface Water Storage','kg/m^2',&
1059	2585a2765,2772
1060	> CALL histdef(hist_id, 'DelSWE', 'Change in Snow Water Equivalent', 'kg/m^2', &
1061	> & iim,jjm, hori_id, 1, 1, 1, -99, 32, sumscatter(1), dt,dw)
1062	> IF ( control_flags%do_irrigation ) THEN
1063	> CALL histdef(hist_id, 'Qirrig', 'Irrigation', 'kg/m^2/s', &
1064	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_scinsec(1), dt,dw)
1065	> CALL histdef(hist_id, 'Qirrig_req', 'Irrigation requirement', 'kg/m^2/s', &
1066	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_scinsec(1), dt,dw)
1067	> ENDIF
1068	2590a2778,2779
1069	> CALL histdef(hist_id, 'PotSurfT', 'Potential (Unstressed) surface temperature', 'K', &
1070	> & iim,jjm, hori_id, 1,1,1, -99, 32, ave(1), dt,dw)
1071	2595c2784,2788
1072	< CALL histdef(hist_id, 'SWE', '3D soil water equivalent','kg/m^2', &
1073	---
1074	> CALL histdef(hist_id, 'SWI', 'Soil wetness index','-', &
1075	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1076	> CALL histdef(hist_id, 'SurfStor', 'Surface Water Storage','kg/m^2', &
1077	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1078	> CALL histdef(hist_id, 'SWE', 'Snow Water Equivalent', 'kg/m^2', &
1079	2615a2809,2810
1080	> CALL histdef(hist_id, 'PotEvapOld', 'Potential evapotranspiration old method', 'kg/m^2/s', &
1081	> & iim,jjm, hori_id, 1, 1, 1, -99, 32, fluxop_scinsec(1), dt,dw)
1082	2621a2817,2818
1083	> CALL histdef(hist_id, 'EWater', 'Open water evaporation', 'kg/m^2/s', &
1084	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_scinsec(1), dt,dw)
1085	2627a2825,2835
1086	> IF ( control_flags%do_floodplains ) THEN
1087	> CALL histdef(hist_id, 'Qflood', 'Floodplain Evaporation', 'kg/m^2/s', &
1088	> & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_scinsec(1), dt,dw)
1089	> ENDIF
1090	> !-
1091	> !- Surface turbulence
1092	> !-
1093	> CALL histdef(hist_id, 'Z0', 'Roughness height', 'm', &
1094	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1095	> CALL histdef(hist_id, 'EffectHeight', 'Effective displacement height (h-d)', 'm', &
1096	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1097	2641,2642c2849,2850
1098	< CALL histdef(hist_id, 'IrrigationMap', 'Map of irrigated areas', 'm^2', &
1099	< & iim,jjm, hori_id, 1,1,1, -99, 32, once(7), dt,dw)
1100	---
1101	> CALL histdef(hist_id, 'SwampMap', 'Map of swamp areas', 'm^2', &
1102	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(1), dt,dw)
1103	2645,2646c2853,2862
1104	< CALL histdef(hist_id, 'humrel', 'Soil moisture stress', '-', &
1105	< & iim,jjm, hori_id, nvm,1,nvm, vegax_id, 32, avescatter(8), dt,dw)
1106	---
1107	> IF ( control_flags%do_irrigation ) THEN
1108	> CALL histdef(hist_id, 'IrrigationMap', 'Map of irrigated areas', 'm^2', &
1109	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(1), dt,dw)
1110	> ENDIF
1111	> IF ( control_flags%do_floodplains ) THEN
1112	> CALL histdef(hist_id, 'FloodplainsMap', 'Map of flooded areas', 'm^2', &
1113	> & iim,jjm, hori_id, 1,1,1, -99, 32, once(1), dt,dw)
1114	> CALL histdef(hist_id, 'FloodFrac', 'Floodplain Fraction', '-', &
1115	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1116	> ENDIF
1117	2663a2880,2881
1118	> CALL histdef (hist_id,'CO2FLUX','Total output CO2 flux', 'gC/day/(m^2 tot)', &
1119	> & iim,jjm, hori_id, 1,1,1, -99,32, ave(1), dt, dw)
1120	2667a2886,2891
1121	> !- Forcing and grid information
1122	> !-
1123	> CALL histdef(hist_id, 'Tair', 'Near surface air temperature at forcing level', 'K', &
1124	> & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
1125	> CALL histdef(hist_id, 'Qair', 'Near surface specific humidity at forcing level', 'g/g', &
1126	> & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
1127	2971a3196,3199
1128	> !
1129	> CALL histdef(hist2_id, 'frac_bare', 'Bare soil fraction for each tile', '-', &
1130	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, avescatter2(3), dt,dw2)
1131	> !
1132	2986a3215,3218
1133	> CALL histdef(hist2_id, 'njsc', 'Soil class used for hydrology', '-', &
1134	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, once2(1), dt,dw2)
1135	> CALL histdef(hist2_id, 'SWI', 'Soil wetness index','-', &
1136	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt,dw2)
1137	3004a3237,3238
1138	> CALL histdef (hist2_id,'CO2FLUX','Total output CO2 flux', 'gC/day/(m^2 tot)', &
1139	> & iim,jjm, hori_id2, 1,1,1, -99,32, fluxop2(5), dt, dw2)
1140	3041a3276,3281
1141	> CALL histdef(hist2_id, 'floodr', 'Floodplains reservoir', 'kg/m^2', &
1142	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(8), dt,dw2)
1143	> CALL histdef(hist2_id, 'floodh', 'Floodplains height', 'mm', &
1144	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(8), dt,dw2)
1145	> CALL histdef(hist2_id, 'pondr', 'Ponds reservoir', 'kg/m^2', &
1146	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(8), dt,dw2)
1147	3045a3286,3287
1148	> CALL histdef(hist2_id, 'returnflow', 'Returnflow to bottom layer', 'mm/d', &
1149	> & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop2(8), dt,dw2)
1150	3049a3292,3295
1151	> CALL histdef(hist2_id, 'floodmap', 'Map of floodplains', 'm^2', &
1152	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(8), dt,dw2)
1153	> CALL histdef(hist2_id, 'swampmap', 'Map of swamps', 'm^2', &
1154	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(8), dt,dw2)
1155	3070c3316,3318
1156	< CALL histdef(hist2_id, 'soiltype', 'Fraction of soil textures', '%', &
1157	---
1158	> CALL histdef(hist2_id, 'vbeta5', 'Beta for bare soil', '-', &
1159	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(9), dt, dw2)
1160	> CALL histdef(hist2_id, 'soiltile', 'Fraction of soil tile', '%', &
1161	3071a3320,3325
1162	> CALL histdef(hist2_id, 'reinf_slope', 'Slope index for each grid box', '-', &
1163	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(9), dt,dw2)
1164	> CALL histdef(hist2_id, 'soilindex', 'Soil index', '-', &
1165	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, once2(9), dt,dw2)
1166	> CALL histdef(hist2_id, 'soildepth', 'Soil Depth', 'm', &
1167	> & iim,jjm, hori_id2, nstm, 1, nstm, soltax_id, 32, once2(1), dt,dw2)
1168	3124c3378
1169	< & nslm, solay, solayax_id2)
1170	---
1171	> & nslm, diaglev(1:nslm), solayax_id2)
1172	3130c3384
1173	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, avescatter2(3), dt, dw2)
1174	---
1175	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, avescatter2(5), dt, dw2)
1176	3132c3386
1177	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, avescatter2(3), dt, dw2)
1178	---
1179	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, avescatter2(5), dt, dw2)
1180	3134c3388
1181	< & iim,jjm, hori_id2, nnobio, 1, nnobio, nobioax_id2, 32, avescatter2(3), dt, dw2)
1182	---
1183	> & iim,jjm, hori_id2, nnobio, 1, nnobio, nobioax_id2, 32, avescatter2(5), dt, dw2)
1184	3153c3407
1185	< & iim,jjm, hori_id2, 1,1,1, -99, 32, ave2(2), dt, dw2)
1186	---
1187	> & iim,jjm, hori_id2, 1,1,1, -99, 32, ave2(5), dt, dw2)
1188	3155c3409
1189	< & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt, dw2)
1190	---
1191	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(5), dt, dw2)
1192	3157c3411
1193	< & iim,jjm, hori_id2, 1,1,1, -99, 32, sumscatter2(1), dt, dw2)
1194	---
1195	> & iim,jjm, hori_id2, 1,1,1, -99, 32, sumscatter2(7), dt, dw2)
1196	3159c3413
1197	< & iim,jjm, hori_id2, 1,1,1, -99, 32, sumscatter2(1), dt, dw2)
1198	---
1199	> & iim,jjm, hori_id2, 1,1,1, -99, 32, sumscatter2(7), dt, dw2)
1200	3176,3178c3430,3432
1201	< & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(1), dt, dw2)
1202	< CALL histdef(hist2_id, 'DelSWE', 'Change in SWE','kg/m^2',&
1203	< & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(1), dt, dw2)
1204	---
1205	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(7), dt, dw2)
1206	> CALL histdef(hist2_id, 'DelSurfStor', 'Change in Surface Water Storage','kg/m^2',&
1207	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(7), dt,dw2)
1208	3180c3434,3436
1209	< & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(1), dt, dw2)
1210	---
1211	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(7), dt, dw2)
1212	> CALL histdef(hist2_id, 'DelSWE', 'Change in interception storage Snow Water Equivalent', 'kg/m^2', &
1213	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, sumscatter2(7), dt, dw2)
1214	3187c3443
1215	< & iim,jjm, hori_id2, 1,1,1, -99, 32, ave2(1), dt, dw2)
1216	---
1217	> & iim,jjm, hori_id2, 1,1,1, -99, 32, ave2(5), dt, dw2)
1218	3188a3445,3446
1219	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(7), dt, dw2)
1220	> CALL histdef(hist2_id, 'SWI', 'Soil wetness index','-', &
1221	3190c3448
1222	< CALL histdef(hist2_id, 'SWE', '3D soil water equivalent','kg/m^2', &
1223	---
1224	> CALL histdef(hist2_id, 'SurfStor', 'Surface Water Storage','kg/m^2', &
1225	3191a3450,3451
1226	> CALL histdef(hist2_id, 'SWE', 'Snow Water Equivalent', 'kg/m^2', &
1227	> & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
1228	3197c3457
1229	< & iim,jjm, hori_id2, 1, 1, 1, -99, 32, avescatter2(1), dt, dw2)
1230	---
1231	> & iim,jjm, hori_id2, 1, 1, 1, -99, 32, avescatter2(7), dt, dw2)
1232	3200c3460
1233	< & iim,jjm, hori_id2, nslm, 1, nslm, solayax_id2, 32, avescatter2(1), dt, dw2)
1234	---
1235	> & iim,jjm, hori_id2, nslm, 1, nslm, solayax_id2, 32, avescatter2(7), dt, dw2)
1236	3205c3465
1237	< & iim,jjm, hori_id2, ngrnd, 1, ngrnd, solax_id2, 32, avescatter2(1), dt, dw2)
1238	---
1239	> & iim,jjm, hori_id2, ngrnd, 1, ngrnd, solax_id2, 32, avescatter2(7), dt, dw2)
1240	3212c3472
1241	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(1), dt, dw2)
1242	---
1243	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(5), dt, dw2)
1244	3214c3474
1245	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(1), dt, dw2)
1246	---
1247	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(5), dt, dw2)
1248	3216c3476,3478
1249	< & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop_scinsec2(1), dt, dw2)
1250	---
1251	> & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop_scinsec2(5), dt, dw2)
1252	> CALL histdef(hist2_id, 'EWater', 'Open water evaporation', 'kg/m^2/s', &
1253	> & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop_scinsec2(5), dt, dw2)
1254	3222c3484
1255	< & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt, dw2)
1256	---
1257	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(5), dt, dw2)
1258	3228c3490
1259	< & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt, dw2)
1260	---
1261	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(7), dt, dw2)
1262	3230c3492
1263	< & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt, dw2)
1264	---
1265	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(7), dt, dw2)
1266	3232c3494
1267	< & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(1), dt, dw2)
1268	---
1269	> & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(7), dt, dw2)
1270	3237c3499,3503
1271	< & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(7), dt, dw2)
1272	---
1273	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(1), dt, dw2)
1274	> CALL histdef(hist2_id, 'FloodplainsMap', 'Map of flooded areas', 'm^2', &
1275	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(1), dt,dw2)
1276	> CALL histdef(hist2_id, 'SwampMap', 'Map of swamp areas', 'm^2', &
1277	> & iim,jjm, hori_id2, 1,1,1, -99, 32, once2(1), dt,dw2)
1278	3239,3241c3505
1279	< & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop_scinsec2(2), dt, dw2)
1280	< CALL histdef(hist2_id, 'humrel', 'Soil moisture stress', '-', &
1281	< & iim,jjm, hori_id2, nvm,1,nvm, vegax_id2, 32, avescatter2(9), dt, dw2)
1282	---
1283	> & iim,jjm, hori_id2, 1,1,1, -99, 32, fluxop_scinsec2(1), dt,dw2)
1284	3247c3511
1285	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(1), dt, dw2)
1286	---
1287	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(3), dt, dw2)
1288	3251c3515
1289	< & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(1), dt, dw2)
1290	---
1291	> & iim,jjm, hori_id2, nvm, 1, nvm, vegax_id2, 32, fluxop_scinsec2(3), dt, dw2)
1292	3258a3523,3525
1293	> ! Total output CO2 flux
1294	> CALL histdef (hist2_id,'CO2FLUX','Total output CO2 flux', 'gC/day/(m^2 tot)', &
1295	> & iim,jjm, hori_id2, 1,1,1, -99,32, ave2(1), dt, dw2)
1296	3368a3636
1297	>
1298	_______________________________________________________________________________________________________________________
1299	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_sechiba/slowproc.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba/slowproc.f90
1300	4c4
1301	< !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/slowproc.f90,v 1.43 2009/08/06 07:52:23 ssipsl Exp $
1302	---
1303	> !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/slowproc.f90,v 1.43.2.1 2009/11/10 17:53:58 ssipsl Exp $
1304	31,41d30
1305	< ! To use OLD or NEW iterpollation schemes :
1306	< INTERFACE slowproc_interlai
1307	< MODULE PROCEDURE slowproc_interlai_OLD, slowproc_interlai_NEW
1308	< END INTERFACE
1309	< INTERFACE slowproc_interpol
1310	< MODULE PROCEDURE slowproc_interpol_OLD, slowproc_interpol_NEW
1311	< END INTERFACE
1312	< INTERFACE slowproc_interpol_g
1313	< MODULE PROCEDURE slowproc_interpol_OLD_g, slowproc_interpol_NEW_g
1314	< END INTERFACE
1315	<
1316	50c39,40
1317	<
1318	---
1319	> REAL(r_std), SAVE :: slope_default = 0.1
1320	> !
1321	56d45
1322	< LOGICAL, SAVE :: old_lai = .FALSE. ! Old Lai Map interpolation
1323	58d46
1324	< LOGICAL, SAVE :: old_veget = .FALSE. ! Old veget Map interpolation
1325	61a50
1326	> REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilclass_default
1327	71c60
1328	< IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
1329	---
1330	> IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
1331	76c65
1332	< lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
1333	---
1334	> lai, frac_bare, height, veget, frac_nobio, njsc, veget_max, totfrac_nobio, qsintmax, &
1335	117a107
1336	> INTEGER(i_std), DIMENSION(kjpindex), INTENT(inout) :: njsc !! Soil type map to be created from the input map
1337	118a109
1338	> REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: frac_bare !! Bare soil fraction for each tile
1339	124c115,116
1340	< REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout):: soiltype !! fraction of soil type
1341	---
1342	> REAL(r_std),DIMENSION (kjpindex,nstm), INTENT(inout):: soiltile !! fraction of soil type
1343	> REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: reinf_slope !! slope coef for reinfiltration
1344	157,158c149,151
1345	< & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,&
1346	< & veget_update, veget_year)
1347	---
1348	> & rest_id, read_lai, lai, frac_bare, veget, frac_nobio, totfrac_nobio, soiltile, reinf_slope, veget_max, njsc, &
1349	> & height, lcanop, veget_update, veget_year)
1350	>
1351	244a238,240
1352	> var_name= 'frac_bare'
1353	> CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, frac_bare, 'scatter', nbp_glo, index_g)
1354	> !
1355	248,249c244,251
1356	< var_name= 'soiltype_frac'
1357	< CALL restput_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, soiltype, 'scatter', nbp_glo, index_g)
1358	---
1359	> var_name= 'soiltile_frac'
1360	> CALL restput_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, soiltile, 'scatter', nbp_glo, index_g)
1361	> !
1362	> var_name= 'njsc'
1363	> CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, REAL(njsc, r_std), 'scatter', nbp_glo, index_g)
1364	> !
1365	> var_name= 'reinf_slope'
1366	> CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, reinf_slope, 'scatter', nbp_glo, index_g)
1367	294c296
1368	< IF (check) WRITE(,) "slowproc: day_counter 3",day_counter
1369	---
1370	> IF (check) WRITE(numout,*) "slowproc: day_counter 3",day_counter
1371	308c310
1372	< IF (check) WRITE(,) "slowproc: day_counter 2",day_counter
1373	---
1374	> IF (check) WRITE(numout,*) "slowproc: day_counter 2",day_counter
1375	360a363,365
1376	> IF (EndOfYear .AND. .NOT. land_use_updated) THEN
1377	> lcchange=.FALSE.
1378	> ENDIF
1379	410c415
1380	< lalo,resolution,lai,month,day,read_lai,laimap)
1381	---
1382	> lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
1383	450,451c455,456
1384	< & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,&
1385	< & veget_update, veget_year)
1386	---
1387	> & rest_id, read_lai, lai, frac_bare, veget, frac_nobio, totfrac_nobio, soiltile, reinf_slope, veget_max, njsc, &
1388	> & height, lcanop, veget_update, veget_year)
1389	472a478
1390	> REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: frac_bare !! Bare soil fraction for each tile
1391	478c484,486
1392	< REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltype !! fraction of soil type
1393	---
1394	> REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! fraction of soil type subvar for hydrological processes
1395	> REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: reinf_slope !! slope coef for reinfiltration
1396	> INTEGER(i_std), DIMENSION(kjpindex), INTENT(out) :: njsc !! Soil type map to be created from the input map
1397	486c494,495
1398	< INTEGER(i_std) :: ji, jv, ier
1399	---
1400	> INTEGER(i_std) :: ji, jv, ier, jst
1401	> !
1402	487a497,498
1403	> LOGICAL :: get_slope
1404	> !
1405	490a502,503
1406	> REAL(r_std), DIMENSION(kjpindex) :: tmp_real
1407	>
1408	491a505
1409	> REAL(r_std), DIMENSION (kjpindex,nscm) :: soilclass !! fraction of soil type
1410	502a517,518
1411	> IF (long_print) WRITE (numout,*) "In slowproc_init"
1412	>
1413	509a526,532
1414	> ALLOCATE (soilclass_default(nscm),stat=ier)
1415	> IF (ier.NE.0) THEN
1416	> WRITE (numout,*) ' error in indexveg allocation. We stop. We need nstm words = ',nscm
1417	> STOP 'hydrol_init'
1418	> END IF
1419	> soilclass_default(:)=undef_sechiba
1420	> !
1421	537a561,563
1422	> !
1423	> IF (long_print) WRITE (numout,*) 'slowproc_init : End of allocations'
1424	>
1425	556a583,584
1426	> IF (long_print) WRITE (numout,*) 'slowproc_init : End of allocationsGot restart file set'
1427	>
1428	612c640
1429	< CALL getin('VEGET_YEAR', veget_year_orig)
1430	---
1431	> CALL getin_p('VEGET_YEAR', veget_year_orig)
1432	655c683,685
1433	< var_name= 'soiltype_frac'
1434	---
1435	> IF (long_print) WRITE (numout,*) 'slowproc_init : End of Land_Use configuration'
1436	> !
1437	> var_name= 'soiltile_frac'
1438	658c688,693
1439	< CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., soiltype, "gather", nbp_glo, index_g)
1440	---
1441	> CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., soiltile, "gather", nbp_glo, index_g)
1442	> !
1443	> var_name= 'reinf_slope'
1444	> CALL ioconf_setatt('UNITS', '-')
1445	> CALL ioconf_setatt('LONG_NAME','Slope coef for reinfiltration')
1446	> CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinf_slope, "gather", nbp_glo, index_g)
1447	664a700,709
1448	> var_name= 'njsc'
1449	> CALL ioconf_setatt('UNITS', '-')
1450	> CALL ioconf_setatt('LONG_NAME','Index of soil type')
1451	> CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tmp_real, "gather", nbp_glo, index_g)
1452	> WHERE ( tmp_real .LT. val_exp )
1453	> njsc = NINT(tmp_real)
1454	> ENDWHERE
1455	> !
1456	> IF (long_print) WRITE (numout,*) 'slowproc_init : End CWRR configuration'
1457	> !
1458	669a715,719
1459	> var_name= 'frac_bare'
1460	> CALL ioconf_setatt('UNITS', '-')
1461	> CALL ioconf_setatt('LONG_NAME','Bare soil fraction')
1462	> CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., frac_bare, "gather", nbp_glo, index_g)
1463	> !
1464	703,718c753,754
1465	<
1466	< !MM, T. d'O. : before in constantes_soil :
1467	< ! diaglev = &
1468	< ! & (/ 0.001, 0.004, 0.01, 0.018, 0.045, &
1469	< ! & 0.092, 0.187, 0.375, 0.750, 1.5, &
1470	< ! & 2.0 /)
1471	< ! Place here because it is used in sechiba and stomate (then in teststomate)
1472	< ! to be in sechiba when teststomate will have disapeared.
1473	< !MM Problem here with dpu which depends on soil type
1474	< DO jv = 1, nbdl-1
1475	< ! first 2.0 is dpu
1476	< ! second 2.0 is average
1477	< diaglev(jv) = dpu_cste/(2*(nbdl-1) -1) ( ( 2(jv-1) -1) + ( 2(jv) -1) ) / 2.0
1478	< ENDDO
1479	< diaglev(nbdl) = dpu_cste
1480	<
1481	---
1482	> !
1483	> !
1484	782c818,826
1485	<
1486	---
1487	> !
1488	> !Config Key = 'PREF_SOIL_VEG'
1489	> !Config Desc = The soil tile number for each vegetation
1490	> !Config Def = 0.1
1491	> !Config Help = Gives the number of the soil tile on which we will
1492	> !Config put each vegetation. This allows to divide the hydrological column
1493	> pref_soil_veg = (/ 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /)
1494	> CALL getin_p('PREF_SOIL_VEG', pref_soil_veg)
1495	> WRITE(numout, *)' PREF_SOIL_VEG, pref_soil_veg = ', pref_soil_veg
1496	853c897
1497	< CALL setvar_p (soiltype, val_exp, 'SOIL_FRACTIONS', soiltype_default)
1498	---
1499	> CALL setvar_p (soilclass, val_exp, 'SOIL_FRACTIONS', soilclass_default)
1500	862a907,915
1501	> !Config Key = REINF_SLOPE
1502	> !Config Desc = Slope coef for reinfiltration
1503	> !Config Def = 0.1
1504	> !Config If = IMPOSE_VEG
1505	> !Config Help = Determines the reinfiltration ratio in the grid box due to flat areas
1506	> !
1507	> slope_default=0.1
1508	> CALL setvar_p (reinf_slope, val_exp, 'SLOPE', slope_default)
1509	>
1510	874a928,942
1511	> !
1512	> soilclass=val_exp
1513	> CALL setvar (soilclass, val_exp, 'SOIL_FRACTIONS', soilclass_default)
1514	> njsc(:) = 0
1515	> soiltile(:,:) = zero
1516	> DO ji = 1, kjpindex
1517	> njsc(ji) = MAXLOC(soilclass(ji,:),1)
1518	> soiltile(:,1) = SUM(frac_nobio(ji,:))
1519	> ENDDO
1520	> DO jv = 1, nvm
1521	> jst = pref_soil_veg(jv)
1522	> DO ji = 1, kjpindex
1523	> soiltile(ji,jst) = soiltile(ji,jst) + veget(ji,jv)
1524	> ENDDO
1525	> ENDDO
1526	884,894d951
1527	<
1528	< !Config Key = SLOWPROC_VEGET_OLD_INTERPOL
1529	< !Config Desc = Flag to use old "interpolation" of vegetation map.
1530	< !Config If = NOT IMPOSE_VEG and NOT LAND_USE
1531	< !Config Def = FALSE
1532	< !Config Help = If you want to recover the old (ie orchidee_1_2 branch)
1533	< !Config "interpolation" of vegetation map.
1534	< !
1535	< old_veget = .FALSE.
1536	< CALL getin_p('SLOWPROC_VEGET_OLD_INTERPOL',old_veget)
1537	<
1538	897,898c954
1539	< IF ( .NOT. old_veget ) THEN
1540	< ! NEW slowproc interpol :
1541	---
1542	> ! slowproc interpol :
1543	900,903d955
1544	< ELSE
1545	< ! OLD slowproc interpol :
1546	< CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, veget_max_g, frac_nobio_g)
1547	< ENDIF
1548	986,996d1037
1549	<
1550	< !Config Key = SLOWPROC_LAI_OLD_INTERPOL
1551	< !Config Desc = Flag to use old "interpolation" of LAI
1552	< !Config If = LAI_MAP
1553	< !Config Def = FALSE
1554	< !Config Help = If you want to recover the old (ie orchidee_1_2 branch)
1555	< !Config "interpolation" of LAI map.
1556	< !
1557	< old_lai = .FALSE.
1558	< CALL getin_p('SLOWPROC_LAI_OLD_INTERPOL',old_lai)
1559	<
1560	1002,1003d1042
1561	< IF ( .NOT. old_lai ) THEN
1562	< ! NEW slowproc interlai :
1563	1005,1008d1043
1564	< ELSE
1565	< ! OLD slowproc interlai :
1566	< CALL slowproc_interlai(kjpindex, lalo, resolution, laimap)
1567	< ENDIF
1568	1014c1049
1569	< lalo,resolution,lai,month,day,read_lai,laimap)
1570	---
1571	> lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
1572	1030,1031d1064
1573	< IF ( .NOT. old_lai ) THEN
1574	< ! NEW slowproc interlai :
1575	1033,1036c1066
1576	< ELSE
1577	< ! OLD slowproc interlai :
1578	< CALL slowproc_interlai(kjpindex, lalo, resolution, laimap)
1579	< ENDIF
1580	---
1581	> !
1582	1042c1072
1583	< lalo,resolution,lai,month,day,read_lai,laimap)
1584	---
1585	> lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
1586	1061c1091,1099
1587	< IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.&
1588	---
1589	> IF ( MINVAL(soilclass) .EQ. MAXVAL(soilclass) .AND. MAXVAL(soilclass) .EQ. val_exp .OR.&
1590	> & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN
1591	>
1592	> CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
1593	>
1594	> ENDIF
1595	>
1596	> IF ( MINVAL(soiltile) .EQ. MAXVAL(soiltile) .AND. MAXVAL(soiltile) .EQ. val_exp .OR.&
1597	> & MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int .OR.&
1598	1062a1101
1599	> CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
1600	1064c1103,1141
1601	< CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction)
1602	---
1603	> ! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
1604	> ! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
1605	> ! in the other modules to perform separated energy balances
1606	> njsc(:) = 0
1607	> soiltile(:,:) = zero
1608	> DO ji = 1, kjpindex
1609	> njsc(ji) = MAXLOC(soilclass(ji,:),1)
1610	> soiltile(:,1) = SUM(frac_nobio(ji,:))
1611	> ENDDO
1612	> DO jv = 1, nvm
1613	> jst = pref_soil_veg(jv)
1614	> DO ji = 1, kjpindex
1615	> soiltile(ji,jst) = soiltile(ji,jst) + veget(ji,jv)
1616	> ENDDO
1617	> ENDDO
1618	>
1619	> ! Avoid soiltile < 0.01
1620	> DO jst = 1, nstm
1621	> DO ji = 1, kjpindex
1622	> IF (soiltile(ji,jst) .LT. 0.01) THEN
1623	> soiltile(ji,MAXLOC(soiltile(ji,:),1)) = soiltile(ji,MAXLOC(soiltile(ji,:),1)) + soiltile(ji,jst)
1624	> soiltile(ji,jst) = zero
1625	> ENDIF
1626	> ENDDO
1627	> ENDDO
1628	> ENDIF
1629	>
1630	>
1631	> !Config Key = GET_SLOPE
1632	> !Config Desc = Read slopes from file and do the interpolation
1633	> !Config Def = .FALSE.
1634	> !Config Help = Needed for reading the slopesfile and doing the interpolation. This will be
1635	> ! used by the re-infiltration parametrization
1636	> get_slope = .FALSE.
1637	> CALL getin_p('GET_SLOPE',get_slope)
1638	>
1639	> IF ( MINVAL(reinf_slope) .EQ. MAXVAL(reinf_slope) .AND. MAXVAL(reinf_slope) .EQ. val_exp .OR. get_slope) THEN
1640	>
1641	> CALL slowproc_slope(kjpindex, lalo, neighbours, resolution, contfrac, reinf_slope)
1642	1069,1096d1145
1643	< !
1644	< ! Select the preferences for the distribution of the 13 PFTs on the 3 soil types.
1645	< !
1646	< vegsoil_dist='EQUI'
1647	< !
1648	< SELECT CASE(vegsoil_dist)
1649	< !
1650	< ! A reasonable choice
1651	< !
1652	< CASE('MAXR')
1653	< pref_soil_veg(:,1) = (/ 1, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 /)
1654	< pref_soil_veg(:,2) = (/ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 /)
1655	< pref_soil_veg(:,3) = (/ 3, 1, 1, 1, 1, 1 ,1 ,1 ,1 ,1 ,1 ,1, 1 /)
1656	< !
1657	< ! Current default : equidistribution.
1658	< !
1659	< CASE('EQUI')
1660	< !
1661	< pref_soil_veg(:,:) = zero
1662	< !
1663	< CASE DEFAULT
1664	< !
1665	< WRITE(numout,*) 'The vegetation soil type preference you have chosen does not exist'
1666	< WRITE(numout,*) 'You chose :', vegsoil_dist
1667	< STOP 'slowproc_init'
1668	< !
1669	< END SELECT
1670	< !
1671	1103a1153,1155
1672	> !
1673	> ! Clean-up
1674	> !
1675	1122a1175,1178
1676	> IF ( ALLOCATED (soilclass_default)) DEALLOCATE (soilclass_default)
1677	> !
1678	>
1679	> !
1680	1489c1545
1681	< SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,mm,dd,read_lai,laimap)
1682	---
1683	> SUBROUTINE slowproc_lai (kjpindex, lcanop, stempdiag, lalo, resolution, lai, frac_bare, mm, dd, read_lai, laimap)
1684	1505a1562
1685	> REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_bare !! Bare fraction for each tile
1686	1508a1566
1687	> REAL(r_std), SAVE :: extcoef !! Extinction coefficient for bare soil with LAI
1688	1509a1568,1577
1689	> !Config Key = EXT_COEF
1690	> !Config Desc = Extinction coefficient for bare soil with LAI
1691	> !Config If =
1692	> !Config Def = 0.5
1693	> !Config Help = will impact frac_bare
1694	> !
1695	> extcoef = 3.
1696	> CALL getin_p('EXT_COEF',extcoef)
1697	> !
1698	> IF ( .NOT. read_lai ) THEN
1699	1511,1515c1579
1700	< ! Test Nathalie. On impose LAI PFT 1 a 0
1701	< ! On boucle sur 2,nvm au lieu de 1,nvm
1702	< lai(: ,1) = 0.0
1703	< DO jv = 2,nvm
1704	< !!$ DO jv = 1,nvm
1705	---
1706	> DO jv = 1,nvm
1707	1523d1586
1708	< IF ( .NOT. read_lai ) THEN
1709	1525,1529d1587
1710	< ELSE
1711	< DO ji = 1,kjpindex
1712	< lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
1713	< ENDDO
1714	< ENDIF
1715	1535d1592
1716	< IF ( .NOT. read_lai ) THEN
1717	1538a1596,1608
1718	> !
1719	> CASE default
1720	> !
1721	> ! 3. Problem
1722	> !
1723	> WRITE (numout,*) 'This kind of lai choice is not possible. '// &
1724	> ' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
1725	> STOP 'slowproc_lai'
1726	>
1727	> END SELECT
1728	>
1729	> ENDDO
1730	> !
1731	1539a1610,1614
1732	> !
1733	> DO jv = 1,nvm
1734	> !
1735	> ! If January
1736	> !
1737	1545a1621,1623
1738	> !
1739	> ! If December
1740	> !
1741	1551a1630,1632
1742	> !
1743	> ! ELSE
1744	> !
1745	1559,1563c1640
1746	< ENDIF
1747	< !
1748	< CASE default
1749	< !
1750	< ! 3. Problem
1751	---
1752	> ENDDO
1753	1565,1569c1642
1754	< WRITE (numout,*) 'This kind of lai choice is not possible. '// &
1755	< ' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
1756	< STOP 'slowproc_lai'
1757	<
1758	< END SELECT
1759	---
1760	> ENDIF
1761	1570a1644,1648
1762	> frac_bare(:,:) = zero
1763	> frac_bare(:,1) = un
1764	> IF (extcoef .LT. 100) THEN
1765	> DO jv=2,nvm
1766	> frac_bare(:,jv) = EXP(-extcoef * lai(:,jv))
1767	1571a1650
1768	> ENDIF
1769	1574a1654
1770	> !!
1771	1576d1655
1772	< !MM TAG 1.6 model !
1773	1578c1657
1774	< SUBROUTINE slowproc_interlai_OLD(nbpt, lalo, resolution, laimap)
1775	---
1776	> SUBROUTINE slowproc_interlai(nbpt, lalo, resolution, neighbours, contfrac, laimap)
1777	1587a1667,1669
1778	> INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
1779	> REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
1780	> !
1781	1594,1595d1675
1782	< REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8
1783	< !
1784	1598,1604c1678,1682
1785	< INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj
1786	< REAL(r_std) :: lev(1), date, dt, coslat, pi
1787	< INTEGER(i_std) :: itau(1)
1788	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu
1789	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask
1790	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful
1791	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: laimaporig
1792	---
1793	> INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
1794	> REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu
1795	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon
1796	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
1797	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
1798	1606,1615c1684,1691
1799	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
1800	< INTEGER, DIMENSION(nbpt) :: n_origlai
1801	< INTEGER, DIMENSION(nbpt) :: n_found
1802	< REAL(r_std), DIMENSION(nbpt,nvm) :: frac_origlai
1803	<
1804	< CHARACTER(LEN=80) :: meter
1805	< REAL(r_std) :: prog, sumf
1806	< LOGICAL :: found
1807	< INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog
1808	< REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
1809	---
1810	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu
1811	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1812	> !
1813	> REAL(r_std) :: totarea, coslat, lmax, lmin, ldelta
1814	> INTEGER(i_std) :: idi, nbvmax, nix, njx
1815	> CHARACTER(LEN=30) :: callsign
1816	> !
1817	> LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
1818	1617c1693
1819	< pi = 4. * ATAN(1.)
1820	---
1821	> INTEGER :: ALLOC_ERR
1822	1621c1697
1823	< !Config If = !LAI_MAP
1824	---
1825	> !Config If = LAI_MAP
1826	1637,1640c1713,1736
1827	< ALLOCATE(lon_lu(iml))
1828	< ALLOCATE(lat_lu(jml))
1829	< ALLOCATE(laimap_lu(iml,jml,nvm,tml))
1830	< ALLOCATE(mask_lu(iml,jml))
1831	---
1832	> ALLOC_ERR=-1
1833	> ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
1834	> IF (ALLOC_ERR/=0) THEN
1835	> WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
1836	> STOP
1837	> ENDIF
1838	> ALLOC_ERR=-1
1839	> ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
1840	> IF (ALLOC_ERR/=0) THEN
1841	> WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
1842	> STOP
1843	> ENDIF
1844	> ALLOC_ERR=-1
1845	> ALLOCATE(laimap_lu(iml,jml,nvm,tml), STAT=ALLOC_ERR)
1846	> IF (ALLOC_ERR/=0) THEN
1847	> WRITE(numout,*) "ERROR IN ALLOCATION of laimap_lu : ",ALLOC_ERR
1848	> STOP
1849	> ENDIF
1850	> ALLOC_ERR=-1
1851	> ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
1852	> IF (ALLOC_ERR/=0) THEN
1853	> WRITE(numout,*) "ERROR IN ALLOCATION of resol_lu : ",ALLOC_ERR
1854	> STOP
1855	> ENDIF
1856	1642d1737
1857	< WRITE(numout,*) 'slowproc_interlai : Reading the LAI file'
1858	1648c1743,1746
1859	< CALL flinget(fid, 'mask', iml, jml, 0, 0, 0, 1, mask_lu)
1860	---
1861	> !
1862	> WHERE (laimap_lu(:,:,:,:) < zero )
1863	> laimap_lu(:,:,:,:) = zero
1864	> ENDWHERE
1865	1655,1657d1752
1866	< CALL bcast(mask_lu)
1867	< !
1868	< WRITE(numout,*) 'slowproc_interlai : ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml)
1869	1658a1754,1765
1870	> ALLOC_ERR=-1
1871	> ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1872	> IF (ALLOC_ERR/=0) THEN
1873	> WRITE(numout,*) "ERROR IN ALLOCATION of lon : ",ALLOC_ERR
1874	> STOP
1875	> ENDIF
1876	> ALLOC_ERR=-1
1877	> ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1878	> IF (ALLOC_ERR/=0) THEN
1879	> WRITE(numout,*) "ERROR IN ALLOCATION of lat : ",ALLOC_ERR
1880	> STOP
1881	> ENDIF
1882	1660,1672c1767,1768
1883	< ijml=iml*jml
1884	< ALLOCATE(lon_ful(ijml))
1885	< ALLOCATE(lat_ful(ijml))
1886	< ALLOCATE(laimaporig(ijml,nvm,tml))
1887	< ALLOCATE(mask(ijml))
1888	<
1889	< DO i=1,iml
1890	< DO j=1,jml
1891	< iki=(j-1)*iml+i
1892	< lon_ful(iki)=lon_lu(i)
1893	< lat_ful(iki)=lat_lu(j)
1894	< laimaporig(iki,:,:)=laimap_lu(i,j,:,:)
1895	< mask(iki)=mask_lu(i,j)
1896	---
1897	> DO ip=1,iml
1898	> lat(ip,:) = lat_lu(:)
1899	1673a1770,1771
1900	> DO jp=1,jml
1901	> lon(:,jp) = lon_lu(:)
1902	1676,1679c1774,1779
1903	< WHERE ( laimaporig(:,:,:) .LT. 0 )
1904	< laimaporig(:,:,:) = 0.
1905	< ENDWHERE
1906	< !
1907	---
1908	> ALLOC_ERR=-1
1909	> ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1910	> IF (ALLOC_ERR/=0) THEN
1911	> WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
1912	> STOP
1913	> ENDIF
1914	1681,1684c1781
1915	< ALLOCATE(lon_up(nbpt))
1916	< ALLOCATE(lon_low(nbpt))
1917	< ALLOCATE(lat_up(nbpt))
1918	< ALLOCATE(lat_low(nbpt))
1919	---
1920	> ! Consider all points a priori
1921	1686c1783
1922	< DO ib =1, nbpt
1923	---
1924	> mask(:,:) = 0
1925	1688,1690c1785,1786
1926	< ! We find the 4 limits of the grid-box. As we transform the resolution of the model
1927	< ! into longitudes and latitudes we do not have the problem of periodicity.
1928	< ! coslat is a help variable here !
1929	---
1930	> DO ip=1,iml
1931	> DO jp=1,jml
1932	1692c1788,1789
1933	< coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
1934	---
1935	> ! Exclude the points where there is never a LAI value. It is probably
1936	> ! an ocean point.
1937	1694,1695c1791,1793
1938	< lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
1939	< lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
1940	---
1941	> IF ( ANY(laimap_lu(ip,jp,:,:) < 20.) ) THEN
1942	> mask(ip,jp) = 1
1943	> ENDIF
1944	1697c1795
1945	< coslat = pi/180. * R_Earth
1946	---
1947	> ! Resolution in longitude
1948	1699,1700c1797,1804
1949	< lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
1950	< lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
1951	---
1952	> coslat = MAX( COS( lat(ip,jp) * pi/180. ), mincos )
1953	> IF ( ip .EQ. 1 ) THEN
1954	> resol_lu(ip,jp,1) = ABS( lon(ip+1,jp) - lon(ip,jp) ) * pi/180. * R_Earth * coslat
1955	> ELSEIF ( ip .EQ. iml ) THEN
1956	> resol_lu(ip,jp,1) = ABS( lon(ip,jp) - lon(ip-1,jp) ) * pi/180. * R_Earth * coslat
1957	> ELSE
1958	> resol_lu(ip,jp,1) = ABS( lon(ip+1,jp) - lon(ip-1,jp) )/2. * pi/180. * R_Earth * coslat
1959	> ENDIF
1960	1701a1806
1961	> ! Resolution in latitude
1962	1702a1808,1814
1963	> IF ( jp .EQ. 1 ) THEN
1964	> resol_lu(ip,jp,2) = ABS( lat(ip,jp) - lat(ip,jp+1) ) * pi/180. * R_Earth
1965	> ELSEIF ( jp .EQ. jml ) THEN
1966	> resol_lu(ip,jp,2) = ABS( lat(ip,jp-1) - lat(ip,jp) ) * pi/180. * R_Earth
1967	> ELSE
1968	> resol_lu(ip,jp,2) = ABS( lat(ip,jp-1) - lat(ip,jp+1) )/2. * pi/180. * R_Earth
1969	> ENDIF
1970	1705,1739c1817
1971	< lon_up = NINT( lon_up * 1E6 ) * 1E-6
1972	< lon_low = NINT( lon_low * 1E6 ) * 1E-6
1973	< lat_up = NINT( lat_up * 1E6 ) * 1E-6
1974	< lat_low = NINT( lat_low * 1E6 ) * 1E-6
1975	< !
1976	< ! Get the limits of the integration domaine so that we can speed up the calculations
1977	< !
1978	< domaine_lon_min = MINVAL(lon_low)
1979	< domaine_lon_max = MAXVAL(lon_up)
1980	< domaine_lat_min = MINVAL(lat_low)
1981	< domaine_lat_max = MAXVAL(lat_up)
1982	< !
1983	< !!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
1984	< !!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
1985	< !
1986	< ! Ensure that the fine grid covers the whole domain
1987	< WHERE ( lon_ful(:) .LT. domaine_lon_min )
1988	< lon_ful(:) = lon_ful(:) + 360.
1989	< ENDWHERE
1990	< !
1991	< WHERE ( lon_ful(:) .GT. domaine_lon_max )
1992	< lon_ful(:) = lon_ful(:) - 360.
1993	< ENDWHERE
1994	< lon_ful = NINT( lon_ful * 1E6 ) * 1E-6
1995	< lat_ful = NINT( lat_ful * 1E6 ) * 1E-6
1996	< !
1997	< WRITE(numout,*) 'Interpolating the lai map :'
1998	< WRITE(numout,'(2a40)')'0%--------------------------------------', &
1999	< & '------------------------------------100%'
2000	< !
2001	< ilast = 1
2002	< n_origlai(:) = 0
2003	< laimap(:,:,:) = 0.
2004	< !
2005	< DO ip=1,ijml
2006	---
2007	> ENDDO
2008	1741c1819,1820
2009	< ! Give a progress meter
2010	---
2011	> ! The number of maximum vegetation map points in the GCM grid is estimated.
2012	> ! Some lmargin is taken.
2013	1743,1751c1822,1825
2014	< ! prog = ip/float(ijml)*79.
2015	< ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(ijml)*79. ) THEN
2016	< ! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
2017	< ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
2018	< ! WRITE(numout, advance="no", FMT='(a80)') meter
2019	< ! ENDIF
2020	< iprog = NINT(float(ip)/float(ijml)79.) - NINT(float(ip-1)/float(ijml)79.)
2021	< IF ( iprog .NE. 0 ) THEN
2022	< WRITE(numout,'(a1,$)') 'y'
2023	---
2024	> IF (is_root_prc) THEN
2025	> nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2
2026	> njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2
2027	> nbvmax = nix*njx
2028	1752a1827
2029	> CALL bcast(nbvmax)
2030	1754,1770c1829
2031	< ! Only start looking for its place in the smaler grid if we are within the domaine
2032	< ! That should speed up things !
2033	< !
2034	< IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
2035	< ( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
2036	< ( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
2037	< ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
2038	< !
2039	< ! look for point on GCM grid which this point on fine grid belongs to.
2040	< ! First look at the point on the model grid where we arrived just before. There is
2041	< ! a good chace that neighbouring points on the fine grid fall into the same model
2042	< ! grid box.
2043	< !
2044	< IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. &
2045	< ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
2046	< ( lat_ful(ip) .GE. lat_low(ilast) ) .AND. &
2047	< ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
2048	---
2049	> callsign = 'LAI map'
2050	1772c1831
2051	< ! We were lucky
2052	---
2053	> ok_interpol = .FALSE.
2054	1774,2052c1833,1834
2055	< IF (mask(ip) .GT. 0) THEN
2056	< n_origlai(ilast) = n_origlai(ilast) + 1
2057	< DO i=1,nvm
2058	< DO j=1,12
2059	< laimap(ilast,i,j) = laimap(ilast,i,j) + laimaporig(ip,i,j)
2060	< ENDDO
2061	< ENDDO
2062	< ENDIF
2063	< !
2064	< ELSE
2065	< !
2066	< ! Otherwise, look everywhere.
2067	< ! Begin close to last grid point.
2068	< !
2069	< found = .FALSE.
2070	< idi = 1
2071	< !
2072	< DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
2073	<
2074	< !
2075	< ! forward and backward
2076	< !
2077	< DO ii = 1,2
2078	< !
2079	< IF ( ii .EQ. 1 ) THEN
2080	< ib = ilast - idi
2081	< ELSE
2082	< ib = ilast + idi
2083	< ENDIF
2084	< !
2085	< IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
2086	< IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. &
2087	< ( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
2088	< ( lat_ful(ip) .GE. lat_low(ib) ) .AND. &
2089	< ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
2090	< !
2091	< IF (mask(ip) .gt. 0) THEN
2092	< DO i=1,nvm
2093	< DO j=1,12
2094	< laimap(ib,i,j) = laimap(ib,i,j) + laimaporig(ip,i,j)
2095	< ENDDO
2096	< ENDDO
2097	< n_origlai(ib) = n_origlai(ib) + 1
2098	< ENDIF
2099	< ilast = ib
2100	< found = .TRUE.
2101	< !
2102	< ENDIF
2103	< ENDIF
2104	< !
2105	< ENDDO
2106	< !
2107	< idi = idi + 1
2108	< !
2109	< ENDDO
2110	< !
2111	< ENDIF ! lucky/not lucky
2112	< !
2113	< ENDIF ! in the domain
2114	< ENDDO
2115	<
2116	<
2117	< ! determine fraction of LAI points in each box of the coarse grid
2118	< DO ip=1,nbpt
2119	< IF ( n_origlai(ip) .GT. 0 ) THEN
2120	< DO jv =1,nvm
2121	< laimap(ip,jv,:) = laimap(ip,jv,:)/REAL(n_origlai(ip),r_std)
2122	< ENDDO
2123	< ELSE
2124	< !
2125	< ! Now we need to handle some exceptions
2126	< !
2127	< IF ( lalo(ip,1) .LT. -56.0) THEN
2128	< ! Antartica
2129	< DO jv =1,nvm
2130	< laimap(ip,jv,:) = 0.
2131	< ENDDO
2132	< !
2133	< ELSE IF ( lalo(ip,1) .GT. 70.0) THEN
2134	< ! Artica
2135	< DO jv =1,nvm
2136	< laimap(ip,jv,:) = 0.
2137	< ENDDO
2138	< !
2139	< ELSE IF ( lalo(ip,1) .GT. 55.0 .AND. lalo(ip,2) .GT. -65.0 .AND. lalo(ip,2) .LT. -20.0) THEN
2140	< ! Greenland
2141	< DO jv =1,nvm
2142	< laimap(ip,jv,:) = 0.
2143	< ENDDO
2144	< !
2145	< ELSE
2146	< !
2147	< WRITE(numout,*) 'PROBLEM, no point in the lai map found for this grid box'
2148	< WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip)
2149	< WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip)
2150	< !
2151	< WRITE(numout,*) 'Looking for nearest point on the lai map file'
2152	< CALL slowproc_nearest (ijml, lon_ful, lat_ful, &
2153	< lalo(ip,2), lalo(ip,1), inear)
2154	< WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', &
2155	< lon_ful(inear),lat_ful(inear)
2156	< !
2157	< DO jv =1,nvm
2158	< laimap(ip,jv,:) = laimaporig(inear,jv,:)
2159	< ENDDO
2160	< ENDIF
2161	< ENDIF
2162	< ENDDO
2163	< !
2164	< WRITE(numout,*) 'slowproc_interlai : Interpolation Done'
2165	< !
2166	< !
2167	< !
2168	< DEALLOCATE(lon_up)
2169	< DEALLOCATE(lon_low)
2170	< DEALLOCATE(lat_up)
2171	< DEALLOCATE(lat_low)
2172	< DEALLOCATE(lat_ful)
2173	< DEALLOCATE(lon_ful)
2174	< DEALLOCATE(lat_lu)
2175	< DEALLOCATE(lon_lu)
2176	< DEALLOCATE(laimap_lu)
2177	< DEALLOCATE(laimaporig)
2178	< DEALLOCATE(mask_lu)
2179	< DEALLOCATE(mask)
2180	< !
2181	< RETURN
2182	< !
2183	< END SUBROUTINE slowproc_interlai_OLD
2184	< !!
2185	< !! Interpolate the LAI map to the grid of the model
2186	< !!
2187	< SUBROUTINE slowproc_interlai_NEW(nbpt, lalo, resolution, neighbours, contfrac, laimap)
2188	< !
2189	< !
2190	< !
2191	< ! 0.1 INPUT
2192	< !
2193	< INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
2194	< REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
2195	< REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
2196	< !
2197	< INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
2198	< REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
2199	< !
2200	< ! 0.2 OUTPUT
2201	< !
2202	< REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned
2203	< !
2204	< ! 0.3 LOCAL
2205	< !
2206	< !
2207	< CHARACTER(LEN=80) :: filename
2208	< INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
2209	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu
2210	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon
2211	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2212	< INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
2213	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu
2214	< INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2215	< !
2216	< REAL(r_std) :: totarea
2217	< INTEGER(i_std) :: idi, nbvmax
2218	< CHARACTER(LEN=30) :: callsign
2219	< !
2220	< LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
2221	< !
2222	< INTEGER :: ALLOC_ERR
2223	< !
2224	< !Config Key = LAI_FILE
2225	< !Config Desc = Name of file from which the vegetation map is to be read
2226	< !Config If = LAI_MAP
2227	< !Config Def = ../surfmap/lai2D.nc
2228	< !Config Help = The name of the file to be opened to read the LAI
2229	< !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2230	< !Config map which is derived from a Nicolas VIOVY one.
2231	< !
2232	< filename = 'lai2D.nc'
2233	< CALL getin_p('LAI_FILE',filename)
2234	< !
2235	< IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
2236	< CALL bcast(iml)
2237	< CALL bcast(jml)
2238	< CALL bcast(lml)
2239	< CALL bcast(tml)
2240	< !
2241	< !
2242	< ALLOC_ERR=-1
2243	< ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
2244	< IF (ALLOC_ERR/=0) THEN
2245	< WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
2246	< STOP
2247	< ENDIF
2248	< ALLOC_ERR=-1
2249	< ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
2250	< IF (ALLOC_ERR/=0) THEN
2251	< WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
2252	< STOP
2253	< ENDIF
2254	< ALLOC_ERR=-1
2255	< ALLOCATE(laimap_lu(iml,jml,nvm,tml), STAT=ALLOC_ERR)
2256	< IF (ALLOC_ERR/=0) THEN
2257	< WRITE(numout,*) "ERROR IN ALLOCATION of laimap_lu : ",ALLOC_ERR
2258	< STOP
2259	< ENDIF
2260	< !
2261	< !
2262	< IF (is_root_prc) THEN
2263	< CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
2264	< CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
2265	< CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu)
2266	< !
2267	< WHERE (laimap_lu(:,:,:,:) < zero )
2268	< laimap_lu(:,:,:,:) = zero
2269	< ENDWHERE
2270	< !
2271	< CALL flinclo(fid)
2272	< ENDIF
2273	< CALL bcast(lon_lu)
2274	< CALL bcast(lat_lu)
2275	< CALL bcast(laimap_lu)
2276	< !
2277	< ALLOC_ERR=-1
2278	< ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
2279	< IF (ALLOC_ERR/=0) THEN
2280	< WRITE(numout,*) "ERROR IN ALLOCATION of lon : ",ALLOC_ERR
2281	< STOP
2282	< ENDIF
2283	< ALLOC_ERR=-1
2284	< ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
2285	< IF (ALLOC_ERR/=0) THEN
2286	< WRITE(numout,*) "ERROR IN ALLOCATION of lat : ",ALLOC_ERR
2287	< STOP
2288	< ENDIF
2289	< !
2290	< DO ip=1,iml
2291	< lat(ip,:) = lat_lu(:)
2292	< ENDDO
2293	< DO jp=1,jml
2294	< lon(:,jp) = lon_lu(:)
2295	< ENDDO
2296	< !
2297	< ALLOC_ERR=-1
2298	< ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2299	< IF (ALLOC_ERR/=0) THEN
2300	< WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
2301	< STOP
2302	< ENDIF
2303	< !
2304	< ! Consider all points a priori
2305	< !
2306	< !!$ mask(:,:) = 1
2307	< mask(:,:) = 0
2308	< !
2309	< DO ip=1,iml
2310	< DO jp=1,jml
2311	< !
2312	< ! Exclude the points where there is never a LAI value. It is probably
2313	< ! an ocean point.
2314	< !
2315	< !!$ IF (ABS(SUM(laimap_lu(ip,jp,:,:))) < EPSILON(laimap_lu) ) THEN
2316	< !!$ mask(ip,jp) = 0
2317	< !!$ ENDIF
2318	< !
2319	< IF ( ANY(laimap_lu(ip,jp,:,:) < 20.) ) THEN
2320	< mask(ip,jp) = 1
2321	< ENDIF
2322	< !
2323	< ENDDO
2324	< ENDDO
2325	< !
2326	< nbvmax = 20
2327	< !
2328	< callsign = 'LAI map'
2329	< !
2330	< ok_interpol = .FALSE.
2331	< DO WHILE ( .NOT. ok_interpol )
2332	< WRITE(numout,*) "Projection arrays for ",callsign," : "
2333	< WRITE(numout,*) "nbvmax = ",nbvmax
2334	---
2335	> WRITE(numout,*) "Projection arrays for ",callsign," : "
2336	> WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
2337	2072d1853
2338	<
2339	2074,2080d1854
2340	< IF ( .NOT. ok_interpol ) THEN
2341	< DEALLOCATE(sub_area)
2342	< DEALLOCATE(sub_index)
2343	< ENDIF
2344	< !
2345	< nbvmax = nbvmax * 2
2346	< ENDDO
2347	2102,2110d1875
2348	<
2349	< !!$ IF ( ANY( laimap(ib,:,:) > 10 ) ) THEN
2350	< !!$ WRITE(numout,*) "For point ",ib
2351	< !!$ WRITE(numout,*) lalo(ib,1), lalo(ib,2)
2352	< !!$ WRITE(numout,*) "For ib=",ib," we have LAI ",laimap(ib,:,1:idi)
2353	< !!$ WRITE(numout,*) "Detail of projection :"
2354	< !!$ WRITE(numout,*) sub_area(ib,1:idi)
2355	< !!$ WRITE(numout,*) "Total for projection :" ,totarea
2356	< !!$ ENDIF
2357	2121a1887,1902
2358	> ! Normelize the read LAI by the values SECHIBA is used to
2359	> !
2360	> DO ib=1,nbpt
2361	> DO jv=1,nvm
2362	> lmax = MAXVAL(laimap(ib,jv,:))
2363	> lmin = MINVAL(laimap(ib,jv,:))
2364	> ldelta = lmax-lmin
2365	> IF ( ldelta < min_sechiba) THEN
2366	> ! LAI constante ... keep it constant
2367	> laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)+(llaimax(jv)+llaimin(jv))/deux
2368	> ELSE
2369	> laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)/(lmax-lmin)*(llaimax(jv)-llaimin(jv))+llaimin(jv)
2370	> ENDIF
2371	> ENDDO
2372	> ENDDO
2373	> !
2374	2133a1915
2375	> DEALLOCATE (resol_lu)
2376	2137c1919
2377	< END SUBROUTINE slowproc_interlai_NEW
2378	---
2379	> END SUBROUTINE slowproc_interlai
2380	2546c2328
2381	< WRITE(,) "SUM(veget_last(",ib,")) = ",sum_veg
2382	---
2383	> WRITE(numout,*) "SUM(veget_last(",ib,")) = ",sum_veg
2384	2718d2499
2385	< !MM TAG 1.6 model !
2386	2720c2501
2387	< SUBROUTINE slowproc_interpol_OLD(nbpt, lalo, neighbours, resolution, veget, frac_nobio )
2388	---
2389	> SUBROUTINE slowproc_interpol(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
2390	2730a2512
2391	> REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
2392	2736a2519,2520
2393	> LOGICAL :: ok_interpol ! optionnal return of aggregate_vec
2394	> !
2395	2739d2522
2396	< REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8
2397	2744,2746c2527
2398	< INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid
2399	< REAL(r_std) :: lev(1), date, dt, coslat, pi
2400	< INTEGER(i_std) :: itau(1)
2401	---
2402	> INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid
2403	2748,2750c2529,2532
2404	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
2405	< INTEGER, DIMENSION(nbpt,nolson) :: n_origveg
2406	< INTEGER, DIMENSION(nbpt) :: n_found
2407	---
2408	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2409	> INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index
2410	> REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg
2411	> REAL(r_std), DIMENSION(nbpt) :: n_found
2412	2754,2758c2536,2540
2413	< CHARACTER(LEN=80) :: meter
2414	< REAL(r_std) :: prog, sumf
2415	< LOGICAL :: found
2416	< INTEGER :: idi, ilast, ii, jv, inear, iprog
2417	< REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
2418	---
2419	> CHARACTER(LEN=40) :: callsign
2420	> REAL(r_std) :: sumf, resol_lon, resol_lat
2421	> INTEGER(i_std) :: idi, jv, inear, nbvmax, nix, njx
2422	> !
2423	> INTEGER :: ALLOC_ERR
2424	2760c2542,2543
2425	< pi = 4. * ATAN(1.)
2426	---
2427	> n_origveg(:,:) = zero
2428	> n_found(:) = zero
2429	2766a2550
2430	> !Config If = !LAND_USE
2431	2783,2792c2567,2591
2432	< ALLOCATE(lat_ful(iml))
2433	< ALLOCATE(lon_ful(iml))
2434	< ALLOCATE(vegmap(iml))
2435	< !
2436	< WRITE(numout,*) 'Reading the vegetation file'
2437	< !
2438	< IF (is_root_prc) THEN
2439	< CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
2440	< CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
2441	< CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
2442	---
2443	> ALLOC_ERR=-1
2444	> ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR)
2445	> IF (ALLOC_ERR/=0) THEN
2446	> WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
2447	> STOP
2448	> ENDIF
2449	> ALLOC_ERR=-1
2450	> ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR)
2451	> IF (ALLOC_ERR/=0) THEN
2452	> WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
2453	> STOP
2454	> ENDIF
2455	> ALLOC_ERR=-1
2456	> ALLOCATE(vegmap(iml), STAT=ALLOC_ERR)
2457	> IF (ALLOC_ERR/=0) THEN
2458	> WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
2459	> STOP
2460	> ENDIF
2461	> !
2462	> WRITE(numout,*) 'Reading the OLSON type vegetation file'
2463	> !
2464	> IF (is_root_prc) THEN
2465	> CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
2466	> CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
2467	> CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
2468	2803,2805c2602,2604
2469	< WRITE(,) 'WARNING -- WARNING'
2470	< WRITE(,) 'The vegetation map has to few vegetation types.'
2471	< WRITE(,) 'If you are lucky it will work but please check'
2472	---
2473	> WRITE(numout,*) 'WARNING -- WARNING'
2474	> WRITE(numout,*) 'The vegetation map has to few vegetation types.'
2475	> WRITE(numout,*) 'If you are lucky it will work but please check'
2476	2807,2808c2606,2607
2477	< WRITE(,) 'More vegetation types in file than the code can'
2478	< WRITE(,) 'deal with.: ', MAXVAL(vegmap), nolson
2479	---
2480	> WRITE(numout,*) 'More vegetation types in file than the code can'
2481	> WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson
2482	2812,2863c2611,2613
2483	< ALLOCATE(lon_up(nbpt))
2484	< ALLOCATE(lon_low(nbpt))
2485	< ALLOCATE(lat_up(nbpt))
2486	< ALLOCATE(lat_low(nbpt))
2487	< !
2488	< DO ib =1, nbpt
2489	< !
2490	< ! We find the 4 limits of the grid-box. As we transform the resolution of the model
2491	< ! into longitudes and latitudes we do not have the problem of periodicity.
2492	< ! coslat is a help variable here !
2493	< !
2494	< coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
2495	< !
2496	< lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
2497	< lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
2498	< !
2499	< coslat = pi/180. * R_Earth
2500	< !
2501	< lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
2502	< lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
2503	< !
2504	< !
2505	< veget(ib,:) = 0.0
2506	< frac_nobio (ib,:) = 0.0
2507	< !
2508	< ENDDO
2509	< !
2510	< ! Get the limits of the integration domaine so that we can speed up the calculations
2511	< !
2512	< domaine_lon_min = MINVAL(lon_low)
2513	< domaine_lon_max = MAXVAL(lon_up)
2514	< domaine_lat_min = MINVAL(lat_low)
2515	< domaine_lat_max = MAXVAL(lat_up)
2516	< !
2517	< !!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
2518	< !!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
2519	< !
2520	< ! Ensure that the fine grid covers the whole domain
2521	< WHERE ( lon_ful(:) .LT. domaine_lon_min )
2522	< lon_ful(:) = lon_ful(:) + 360.
2523	< ENDWHERE
2524	< !
2525	< WHERE ( lon_ful(:) .GT. domaine_lon_max )
2526	< lon_ful(:) = lon_ful(:) - 360.
2527	< ENDWHERE
2528	< !
2529	< WRITE(numout,*) 'Interpolating the vegetation map :'
2530	< WRITE(numout,'(2a40)')'0%--------------------------------------', &
2531	< & '------------------------------------100%'
2532	< !
2533	< ilast = 1
2534	< n_origveg(:,:) = 0
2535	---
2536	> ! Some assumptions on the vegetation file. This information should be
2537	> ! be computed or read from the file.
2538	> ! It is the reolution in meters of the grid of the vegetation file.
2539	2865c2615,2616
2540	< DO ip=1,iml
2541	---
2542	> resol_lon = 5000.
2543	> resol_lat = 5000.
2544	2867c2618,2619
2545	< ! Give a progress meter
2546	---
2547	> ! The number of maximum vegetation map points in the GCM grid is estimated.
2548	> ! Some lmargin is taken.
2549	2869,2877c2621,2624
2550	< ! prog = ip/float(iml)*79.
2551	< ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN
2552	< ! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
2553	< ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
2554	< ! WRITE(numout, advance="no", FMT='(a80)') meter
2555	< ! ENDIF
2556	< iprog = NINT(float(ip)/float(iml)79.) - NINT(float(ip-1)/float(iml)79.)
2557	< IF ( iprog .NE. 0 ) THEN
2558	< WRITE(numout,'(a1,$)') 'x'
2559	---
2560	> IF (is_root_prc) THEN
2561	> nix=INT(MAXVAL(resolution_g(:,1))*2/resol_lon)+2
2562	> njx=INT(MAXVAL(resolution_g(:,2))*2/resol_lon)+2
2563	> nbvmax = nix*njx
2564	2878a2626
2565	> CALL bcast(nbvmax)
2566	2880,2916c2628
2567	< ! Only start looking for its place in the smaler grid if we are within the domaine
2568	< ! That should speed up things !
2569	< !
2570	< IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
2571	< ( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
2572	< ( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
2573	< ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
2574	< !
2575	< ! look for point on GCM grid which this point on fine grid belongs to.
2576	< ! First look at the point on the model grid where we arrived just before. There is
2577	< ! a good chace that neighbouring points on the fine grid fall into the same model
2578	< ! grid box.
2579	< !
2580	< !
2581	< ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE
2582	< ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED.
2583	< !
2584	< IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. &
2585	< ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
2586	< ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. &
2587	< ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
2588	< !
2589	< ! We were lucky
2590	< !
2591	< n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1
2592	< !
2593	< ELSE
2594	< !
2595	< ! Otherwise, look everywhere.
2596	< ! Begin close to last grid point.
2597	< !
2598	< found = .FALSE.
2599	< idi = 1
2600	< !
2601	< DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
2602	< !
2603	< ! forward and backward
2604	---
2605	> callsign="Vegetation map"
2606	2918c2630,2633
2607	< DO ii = 1,2
2608	---
2609	> ok_interpol = .FALSE.
2610	>
2611	> WRITE(numout,*) "Projection arrays for ",callsign," : "
2612	> WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
2613	2920,2923c2635,2639
2614	< IF ( ii .EQ. 1 ) THEN
2615	< ib = ilast - idi
2616	< ELSE
2617	< ib = ilast + idi
2618	---
2619	> ALLOC_ERR=-1
2620	> ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
2621	> IF (ALLOC_ERR/=0) THEN
2622	> WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
2623	> STOP
2624	2924a2641,2648
2625	> sub_index(:,:)=0
2626	> ALLOC_ERR=-1
2627	> ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
2628	> IF (ALLOC_ERR/=0) THEN
2629	> WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
2630	> STOP
2631	> ENDIF
2632	> sub_area(:,:)=zero
2633	2926,2934c2650,2651
2634	< IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
2635	< IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. &
2636	< ( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
2637	< ( lat_ful(ip) .GT. lat_low(ib) ) .AND. &
2638	< ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
2639	< !
2640	< n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1
2641	< ilast = ib
2642	< found = .TRUE.
2643	---
2644	> WRITE(numout,*) 'Carteveg range LON:', MINVAL(lon_ful), MAXVAL(lon_ful)
2645	> WRITE(numout,*) 'Carteveg range LAT:', MINVAL(lat_ful), MAXVAL(lat_ful)
2646	2936,2937c2653,2655
2647	< ENDIF
2648	< ENDIF
2649	---
2650	> CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, &
2651	> & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
2652	> & nbvmax, sub_index, sub_area, ok_interpol)
2653	2939d2656
2654	< ENDDO
2655	2940a2658,2663
2656	> DO ib = 1, nbpt
2657	> idi=1
2658	> DO WHILE ( sub_area(ib,idi) > zero )
2659	> ip = sub_index(ib,idi)
2660	> n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi)
2661	> n_found(ib) = n_found(ib) + sub_area(ib,idi)
2662	2942d2664
2663	< !
2664	2944,2947d2665
2665	< !
2666	< ENDIF ! lucky/not lucky
2667	< !
2668	< ENDIF ! in the domain
2669	2949c2667
2670	<
2671	---
2672	> !
2673	2954,2962d2671
2674	<
2675	< !
2676	< ! determine number of points of the fine grid which fall into each box of the
2677	< ! coarse grid
2678	< !
2679	< DO ib = 1, nbpt
2680	< n_found(ib) = SUM( n_origveg(ib,:) )
2681	< ENDDO
2682	<
2683	2968c2677
2684	< frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std)
2685	---
2686	> frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:)
2687	2973d2681
2688	<
2689	2977a2686,2688
2690	> veget(:,:) = zero
2691	> frac_nobio(:,:) = zero
2692	> !
2693	2990,2992c2701
2694	< !
2695	< WRITE(numout,*)
2696	< WRITE(numout,*) 'Interpolation Done'
2697	---
2698	> WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
2699	3024,3026c2733,2735
2700	< WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box'
2701	< WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib)
2702	< WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib)
2703	---
2704	> WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib
2705	> WRITE(numout,*) 'Longitude range : ', lalo(ib,2)
2706	> WRITE(numout,*) 'Latitude range : ', lalo(ib,1)
2707	3062,3067d2770
2708	< DEALLOCATE(lon_up)
2709	< DEALLOCATE(lon_low)
2710	< DEALLOCATE(lat_up)
2711	< DEALLOCATE(lat_low)
2712	< DEALLOCATE(lat_ful)
2713	< DEALLOCATE(lon_ful)
2714	3068a2772,2775
2715	> DEALLOCATE(lat_ful, lon_ful)
2716	> DEALLOCATE(sub_index)
2717	> DEALLOCATE(sub_area)
2718	>
2719	3072c2779
2720	< END SUBROUTINE slowproc_interpol_OLD
2721	---
2722	> END SUBROUTINE slowproc_interpol
2723	3076c2783
2724	< SUBROUTINE slowproc_interpol_NEW(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
2725	---
2726	> SUBROUTINE slowproc_interpol_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
2727	3113c2820
2728	< INTEGER(i_std) :: idi, jv, inear, nbvmax
2729	---
2730	> INTEGER(i_std) :: idi, jv, inear, nbvmax, nix, njx
2731	3133c2840
2732	< CALL getin_p('VEGETATION_FILE',filename)
2733	---
2734	> CALL getin('VEGETATION_FILE',filename)
2735	3135,3139c2842
2736	< if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
2737	< CALL bcast(iml)
2738	< CALL bcast(jml)
2739	< CALL bcast(lml)
2740	< CALL bcast(tml)
2741	---
2742	> CALL flininfo(filename, iml, jml, lml, tml, fid)
2743	3163d2865
2744	< IF (is_root_prc) THEN
2745	3167a2870,2872
2746	> WRITE(numout,*) 'File name : ', filename
2747	> WRITE(numout,*) 'Min and max vegetation numbers : ', MINVAL(vegmap), MAXVAL(vegmap)
2748	> !
2749	3169,3174d2873
2750	< ENDIF
2751	<
2752	< CALL bcast(lon_ful)
2753	< CALL bcast(lat_ful)
2754	< CALL bcast(vegmap)
2755	<
2756	3193,3195c2892,2900
2757	< ! The number of maximum vegetation map points in the GCM grid should
2758	< ! also be computed and not imposed here.
2759	< nbvmax = iml/nbpt
2760	---
2761	> !
2762	> ! The number of maximum vegetation map points in the GCM grid is estimated.
2763	> ! Some margin is taken.
2764	> !
2765	> nix=INT(MAXVAL(resolution_g(:,1)*2)/resol_lon)+1
2766	> njx=INT(MAXVAL(resolution_g(:,2)*2)/resol_lon)+1
2767	> nbvmax = nix*njx
2768	> !
2769	> ! No need to broadcast as this routine is only called on root_proc
2770	3198a2904,2905
2771	> ! We are on a mono proc routine and thus nbvmax does not need to be broadcasted
2772	> !
2773	3200d2906
2774	< DO WHILE ( .NOT. ok_interpol )
2775	3202c2908
2776	< WRITE(numout,*) "nbvmax = ",nbvmax
2777	---
2778	> WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
2779	3219c2925
2780	< CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, &
2781	---
2782	> CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, &
2783	3223,3228d2928
2784	< IF ( .NOT. ok_interpol ) THEN
2785	< DEALLOCATE(sub_area)
2786	< DEALLOCATE(sub_index)
2787	< !
2788	< nbvmax = nbvmax * 2
2789	< ELSE
2790	3240,3241d2939
2791	< ENDIF
2792	< ENDDO
2793	3353c3051,3052
2794	< END SUBROUTINE slowproc_interpol_NEW
2795	---
2796	> END SUBROUTINE slowproc_interpol_g
2797	>
2798	3356,3357c3055,3106
2799	< !! Interpolate the IGBP vegetation map to the grid of the model
2800	< !MM TAG 1.6 model !
2801	---
2802	> !! looks for nearest grid point on the fine map
2803	> !!
2804	> SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
2805	>
2806	> INTEGER(i_std), INTENT(in) :: iml
2807	> REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5
2808	> REAL(r_std), INTENT(in) :: lonmod, latmod
2809	>
2810	> INTEGER(i_std), INTENT(out) :: inear
2811	>
2812	> REAL(r_std) :: pi
2813	> REAL(r_std) :: pa, p
2814	> REAL(r_std) :: coscolat, sincolat
2815	> REAL(r_std) :: cospa, sinpa
2816	> REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
2817	> INTEGER(i_std) :: i
2818	> INTEGER(i_std), DIMENSION(1) :: ineartab
2819	> INTEGER :: ALLOC_ERR
2820	>
2821	> ALLOC_ERR=-1
2822	> ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
2823	> IF (ALLOC_ERR/=0) THEN
2824	> WRITE(numout,*) "ERROR IN ALLOCATION of cosang : ",ALLOC_ERR
2825	> STOP
2826	> ENDIF
2827	>
2828	> pi = 4.0 * ATAN(1.0)
2829	>
2830	> pa = pi/2.0 - latmod*pi/180.0 ! dist. entre pole n et point a
2831	> cospa = COS(pa)
2832	> sinpa = SIN(pa)
2833	>
2834	> DO i = 1, iml
2835	>
2836	> sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 )
2837	> coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 )
2838	>
2839	> p = (lonmod-lon5(i))*pi/180.0 ! angle entre a et b (leurs meridiens)
2840	>
2841	> ! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
2842	> cosang(i) = cospacoscolat + sinpasincolat*COS(p)
2843	>
2844	> ENDDO
2845	>
2846	> ineartab = MAXLOC( cosang(:) )
2847	> inear = ineartab(1)
2848	>
2849	> DEALLOCATE(cosang)
2850	> END SUBROUTINE slowproc_nearest
2851	>
2852	> !!
2853	> !! Interpolate the Zobler soil type map
2854	3359c3108
2855	< SUBROUTINE slowproc_interpol_OLD_g(nbpt, lalo, neighbours, resolution, veget, frac_nobio )
2856	---
2857	> SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
2858	3361a3111,3115
2859	> ! This subroutine should read the Zobler map and interpolate to the model grid. The method
2860	> ! is to get fraction of the three main soiltypes for each grid box.
2861	> ! The soil fraction are going to be put into the array soiltype in the following order :
2862	> ! coarse, medium and fine.
2863	> !
2864	3369a3124
2865	> REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
2866	3373,3374c3128,3129
2867	< REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
2868	< REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
2869	---
2870	> REAL(r_std), INTENT(out) :: soilclass(nbpt, nstm) ! Soil type map to be created from the Zobler map
2871	> REAL(r_std), INTENT(out) :: clayfraction(nbpt) ! The fraction of clay as used by STOMATE
2872	3376d3130
2873	< ! 0.3 LOCAL
2874	3378,3379c3132
2875	< REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8
2876	< INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
2877	---
2878	> ! 0.3 LOCAL
2879	3380a3134
2880	> INTEGER(i_std) :: nbvmax
2881	3383,3384c3137,3138
2882	< INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid
2883	< REAL(r_std) :: lev(1), date, dt, coslat, pi
2884	---
2885	> INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, nbexp
2886	> REAL(r_std) :: lev(1), date, dt
2887	3386,3397c3140,3150
2888	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
2889	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
2890	< INTEGER, DIMENSION(nbpt,nolson) :: n_origveg
2891	< INTEGER, DIMENSION(nbpt) :: n_found
2892	< REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
2893	< REAL(r_std) :: vegcorr(nolson,nvm)
2894	< REAL(r_std) :: nobiocorr(nolson,nnobio)
2895	< CHARACTER(LEN=80) :: meter
2896	< REAL(r_std) :: prog, sumf
2897	< LOGICAL :: found
2898	< INTEGER :: idi, ilast, ii, jv, inear, iprog
2899	< REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
2900	---
2901	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel
2902	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soiltext, soiltext2
2903	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2904	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2905	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
2906	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu
2907	> INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt, solt2
2908	> REAL(r_std) :: sgn, coslat
2909	> CHARACTER(LEN=30) :: callsign
2910	> CHARACTER(LEN=30) :: classif
2911	> INTEGER(i_std) :: nix, njx
2912	3399c3152
2913	< pi = 4. * ATAN(1.)
2914	---
2915	> ! Number of texture classes in Zobler
2916	3401c3154,3155
2917	< CALL get_vegcorr (nolson,vegcorr,nobiocorr)
2918	---
2919	> INTEGER(i_std), PARAMETER :: nzobler = 7
2920	> REAL(r_std),ALLOCATABLE :: textfrac_table(:,:)
2921	3403,3410c3157
2922	< !Config Key = VEGETATION_FILE
2923	< !Config Desc = Name of file from which the vegetation map is to be read
2924	< !Config If = !IMPOSE_VEG
2925	< !Config Def = ../surfmap/carteveg5km.nc
2926	< !Config Help = The name of the file to be opened to read the vegetation
2927	< !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2928	< !Config map which is derived from the IGBP one. We assume that we have
2929	< !Config a classification in 87 types. This is Olson modified by Viovy.
2930	---
2931	> LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
2932	3412,3413c3159
2933	< filename = '../surfmap/carteveg5km.nc'
2934	< CALL getin('VEGETATION_FILE',filename)
2935	---
2936	> INTEGER :: ALLOC_ERR
2937	3415d3160
2938	< CALL flininfo(filename, iml, jml, lml, tml, fid)
2939	3416a3162
2940	> ! Needs to be a configurable variable
2941	3418,3420d3163
2942	< ALLOCATE(lat_ful(iml))
2943	< ALLOCATE(lon_ful(iml))
2944	< ALLOCATE(vegmap(iml))
2945	3422c3165,3173
2946	< WRITE(numout,*) 'Reading the vegetation file'
2947	---
2948	> !Config Key = SOILCLASS_FILE
2949	> !Config Desc = Name of file from which soil types are read
2950	> !Config Def = ../surfmap/soils_param.nc
2951	> !Config If = !IMPOSE_VEG
2952	> !Config Help = The name of the file to be opened to read the soil types.
2953	> !Config The data from this file is then interpolated to the grid of
2954	> !Config of the model. The aim is to get fractions for sand loam and
2955	> !Config clay in each grid box. This information is used for soil hydrology
2956	> !Config and respiration.
2957	3424,3426c3175,3176
2958	< CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
2959	< CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
2960	< CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
2961	---
2962	> filename = '../surfmap/soils_param.nc'
2963	> CALL getin_p('SOILCLASS_FILE',filename)
2964	3428,3429c3178,3184
2965	< CALL flinclo(fid)
2966	<
2967	---
2968	> !Config Key = SOIL_CLASSIF
2969	> !Config Desc = Type of classification used for the map of soil types
2970	> !Config Def = zobler
2971	> !Config If = !IMPOSE_VEG
2972	> !Config Help = The classification used in the file that we use here
2973	> !Config = There are three classification supported:
2974	> !Config = FAO (3 soil types), Zobler (7) and USDA (12)
2975	3431,3438c3186,3192
2976	< IF (MAXVAL(vegmap) .LT. nolson) THEN
2977	< WRITE(,) 'WARNING -- WARNING'
2978	< WRITE(,) 'The vegetation map has to few vegetation types.'
2979	< WRITE(,) 'If you are lucky it will work but please check'
2980	< ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
2981	< WRITE(,) 'More vegetation types in file than the code can'
2982	< WRITE(,) 'deal with.: ', MAXVAL(vegmap), nolson
2983	< STOP 'slowproc_interpol'
2984	---
2985	> classif = 'zobler'
2986	> CALL getin_p('SOIL_CLASSIF',classif)
2987	> !
2988	> !
2989	> IF (is_root_prc) THEN
2990	> CALL flininfo(filename,iml, jml, lml, tml, fid)
2991	> CALL flinclo(fid)
2992	3439a3194,3197
2993	> CALL bcast(iml)
2994	> CALL bcast(jml)
2995	> CALL bcast(lml)
2996	> CALL bcast(tml)
2997	3441,3444c3199
2998	< ALLOCATE(lon_up(nbpt))
2999	< ALLOCATE(lon_low(nbpt))
3000	< ALLOCATE(lat_up(nbpt))
3001	< ALLOCATE(lat_low(nbpt))
3002	---
3003	> ! soils_param.nc file is 1° soit texture file.
3004	3446c3201,3236
3005	< DO ib =1, nbpt
3006	---
3007	> ALLOC_ERR=-1
3008	> ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
3009	> IF (ALLOC_ERR/=0) THEN
3010	> WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
3011	> STOP
3012	> ENDIF
3013	> ALLOC_ERR=-1
3014	> ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
3015	> IF (ALLOC_ERR/=0) THEN
3016	> WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
3017	> STOP
3018	> ENDIF
3019	> ALLOC_ERR=-1
3020	> ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
3021	> IF (ALLOC_ERR/=0) THEN
3022	> WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
3023	> STOP
3024	> ENDIF
3025	> ALLOC_ERR=-1
3026	> ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR)
3027	> IF (ALLOC_ERR/=0) THEN
3028	> WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
3029	> STOP
3030	> ENDIF
3031	> ALLOC_ERR=-1
3032	> ALLOCATE(soiltext2(iml,jml), STAT=ALLOC_ERR)
3033	> IF (ALLOC_ERR/=0) THEN
3034	> WRITE(numout,*) "ERROR IN ALLOCATION of soiltext2 : ",ALLOC_ERR
3035	> STOP
3036	> ENDIF
3037	> ALLOC_ERR=-1
3038	> ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
3039	> IF (ALLOC_ERR/=0) THEN
3040	> WRITE(numout,*) "ERROR IN ALLOCATION of resol_lu : ",ALLOC_ERR
3041	> STOP
3042	> ENDIF
3043	3448,3450c3238,3244
3044	< ! We find the 4 limits of the grid-box. As we transform the resolution of the model
3045	< ! into longitudes and latitudes we do not have the problem of periodicity.
3046	< ! coslat is a help variable here !
3047	---
3048	> IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
3049	> CALL bcast(lon_rel)
3050	> CALL bcast(lat_rel)
3051	> CALL bcast(itau)
3052	> CALL bcast(date)
3053	> CALL bcast(dt)
3054	>
3055	3452c3246,3247
3056	< coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
3057	---
3058	> IF (is_root_prc) CALL flinget(fid, 'soiltext', iml, jml, lml, tml, 1, 1, soiltext)
3059	> CALL bcast(soiltext)
3060	3454,3455c3249,3252
3061	< lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
3062	< lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
3063	---
3064	> IF (classif .EQ. "fao2") THEN
3065	> IF (is_root_prc) CALL flinget(fid, 'soiltext2', iml, jml, lml, tml, 1, 1, soiltext2)
3066	> CALL bcast(soiltext2)
3067	> ENDIF
3068	3457c3254
3069	< coslat = pi/180. * R_Earth
3070	---
3071	> IF (is_root_prc) CALL flinclo(fid)
3072	3459,3460c3256
3073	< lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
3074	< lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
3075	---
3076	> nbexp = 0
3077	3463,3464c3259
3078	< veget(ib,:) = 0.0
3079	< frac_nobio (ib,:) = 0.0
3080	---
3081	> ! Mask of permitted variables.
3082	3466c3261,3263
3083	< ENDDO
3084	---
3085	> mask(:,:) = zero
3086	> DO ip=1,iml
3087	> DO jp=1,jml
3088	3468c3265,3267
3089	< ! Get the limits of the integration domaine so that we can speed up the calculations
3090	---
3091	> IF (soiltext(ip,jp) .GT. min_sechiba) THEN
3092	> mask(ip,jp) = un
3093	> ENDIF
3094	3470,3481c3269
3095	< domaine_lon_min = MINVAL(lon_low)
3096	< domaine_lon_max = MAXVAL(lon_up)
3097	< domaine_lat_min = MINVAL(lat_low)
3098	< domaine_lat_max = MAXVAL(lat_up)
3099	< !
3100	< !!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
3101	< !!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
3102	< !
3103	< ! Ensure that the fine grid covers the whole domain
3104	< WHERE ( lon_ful(:) .LT. domaine_lon_min )
3105	< lon_ful(:) = lon_ful(:) + 360.
3106	< ENDWHERE
3107	---
3108	> ! Resolution in longitude
3109	3483,3485c3271,3278
3110	< WHERE ( lon_ful(:) .GT. domaine_lon_max )
3111	< lon_ful(:) = lon_ful(:) - 360.
3112	< ENDWHERE
3113	---
3114	> coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )
3115	> IF ( ip .EQ. 1 ) THEN
3116	> resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip,jp) ) * pi/180. * R_Earth * coslat
3117	> ELSEIF ( ip .EQ. iml ) THEN
3118	> resol_lu(ip,jp,1) = ABS( lon_rel(ip,jp) - lon_rel(ip-1,jp) ) * pi/180. * R_Earth * coslat
3119	> ELSE
3120	> resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip-1,jp) )/2. * pi/180. * R_Earth * coslat
3121	> ENDIF
3122	3487,3489c3280
3123	< WRITE(numout,*) 'Interpolating the vegetation map :'
3124	< WRITE(numout,'(2a40)')'0%--------------------------------------', &
3125	< & '------------------------------------100%'
3126	---
3127	> ! Resolution in latitude
3128	3491,3492c3282,3288
3129	< ilast = 1
3130	< n_origveg(:,:) = 0
3131	---
3132	> IF ( jp .EQ. 1 ) THEN
3133	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp) - lat_rel(ip,jp+1) ) * pi/180. * R_Earth
3134	> ELSEIF ( jp .EQ. jml ) THEN
3135	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp-1) - lat_rel(ip,jp) ) * pi/180. * R_Earth
3136	> ELSE
3137	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp-1) - lat_rel(ip,jp+1) )/2. * pi/180. * R_Earth
3138	> ENDIF
3139	3494c3290,3291
3140	< DO ip=1,iml
3141	---
3142	> ENDDO
3143	> ENDDO
3144	3496c3293,3294
3145	< ! Give a progress meter
3146	---
3147	> ! The number of maximum vegetation map points in the GCM grid is estimated.
3148	> ! Some lmargin is taken.
3149	3498,3506c3296,3299
3150	< ! prog = ip/float(iml)*79.
3151	< ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN
3152	< ! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
3153	< ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
3154	< ! WRITE(numout, advance="no", FMT='(a80)') meter
3155	< ! ENDIF
3156	< iprog = NINT(float(ip)/float(iml)79.) - NINT(float(ip-1)/float(iml)79.)
3157	< IF ( iprog .NE. 0 ) THEN
3158	< WRITE(numout,'(a1,$)') 'x'
3159	---
3160	> IF (is_root_prc) THEN
3161	> nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2
3162	> njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2
3163	> nbvmax = nix*njx
3164	3507a3301
3165	> CALL bcast(nbvmax)
3166	3509,3510c3303
3167	< ! Only start looking for its place in the smaler grid if we are within the domaine
3168	< ! That should speed up things !
3169	---
3170	> callsign = "Soil types"
3171	3512,3515c3305
3172	< IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
3173	< ( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
3174	< ( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
3175	< ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
3176	---
3177	> ok_interpol = .FALSE.
3178	3517,3520c3307,3308
3179	< ! look for point on GCM grid which this point on fine grid belongs to.
3180	< ! First look at the point on the model grid where we arrived just before. There is
3181	< ! a good chace that neighbouring points on the fine grid fall into the same model
3182	< ! grid box.
3183	---
3184	> WRITE(numout,*) "Projection arrays for ",callsign," : "
3185	> WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
3186	3521a3310,3335
3187	> ALLOC_ERR=-1
3188	> ALLOCATE(solt(nbvmax), STAT=ALLOC_ERR)
3189	> IF (ALLOC_ERR/=0) THEN
3190	> WRITE(numout,*) "ERROR IN ALLOCATION of solt : ",ALLOC_ERR
3191	> STOP
3192	> ENDIF
3193	> ALLOC_ERR=-1
3194	> ALLOCATE(solt2(nbvmax), STAT=ALLOC_ERR)
3195	> IF (ALLOC_ERR/=0) THEN
3196	> WRITE(numout,*) "ERROR IN ALLOCATION of solt2 : ",ALLOC_ERR
3197	> STOP
3198	> ENDIF
3199	> ALLOC_ERR=-1
3200	> ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
3201	> IF (ALLOC_ERR/=0) THEN
3202	> WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3203	> STOP
3204	> ENDIF
3205	> sub_index(:,:,:)=0
3206	> ALLOC_ERR=-1
3207	> ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
3208	> IF (ALLOC_ERR/=0) THEN
3209	> WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3210	> STOP
3211	> ENDIF
3212	> sub_area(:,:)=zero
3213	3523,3524c3337,3339
3214	< ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE
3215	< ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED.
3216	---
3217	> CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
3218	> & iml, jml, lon_rel, lat_rel, mask, callsign, &
3219	> & nbvmax, sub_index, sub_area, ok_interpol)
3220	3526,3529d3340
3221	< IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. &
3222	< ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
3223	< ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. &
3224	< ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
3225	3531c3342
3226	< ! We were lucky
3227	---
3228	> ! Depending on the selected classification, the interpolation will be different
3229	3533c3344,3351
3230	< n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1
3231	---
3232	> SELECTCASE(classif)
3233	> CASE('none')
3234	> ALLOCATE(textfrac_table(nscm,ntext))
3235	> DO ib=1, nbpt
3236	> soilclass(ib,:) = soilclass_default_fao
3237	> clayfraction(ib) = clayfraction_default
3238	> ENDDO
3239	> CASE('zobler')
3240	3535c3353
3241	< ELSE
3242	---
3243	> soilclass_default=soilclass_default_fao
3244	3537,3538c3355
3245	< ! Otherwise, look everywhere.
3246	< ! Begin close to last grid point.
3247	---
3248	> PRINT *, "Using a soilclass map with Zobler classification"
3249	3540,3541c3357
3250	< found = .FALSE.
3251	< idi = 1
3252	---
3253	> ALLOCATE(textfrac_table(nzobler,ntext))
3254	3543c3359
3255	< DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
3256	---
3257	> CALL get_soilcorr (nzobler, textfrac_table)
3258	3545d3360
3259	< ! forward and backward
3260	3547c3362
3261	< DO ii = 1,2
3262	---
3263	> DO ib =1, nbpt
3264	3549,3553c3364
3265	< IF ( ii .EQ. 1 ) THEN
3266	< ib = ilast - idi
3267	< ELSE
3268	< ib = ilast + idi
3269	< ENDIF
3270	---
3271	> ! GO through the point we have found
3272	3555,3559d3365
3273	< IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
3274	< IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. &
3275	< ( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
3276	< ( lat_ful(ip) .GT. lat_low(ib) ) .AND. &
3277	< ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
3278	3561,3563c3367
3279	< n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1
3280	< ilast = ib
3281	< found = .TRUE.
3282	---
3283	> fopt = COUNT(sub_area(ib,:) > zero)
3284	3565,3566c3369
3285	< ENDIF
3286	< ENDIF
3287	---
3288	> ! Check that we found some points
3289	3568c3371,3372
3290	< ENDDO
3291	---
3292	> soilclass(ib,:) = 0.0
3293	> clayfraction(ib) = 0.0
3294	3570c3374,3378
3295	< idi = idi + 1
3296	---
3297	> IF ( fopt .EQ. 0) THEN
3298	> nbexp = nbexp + 1
3299	> soilclass(ib,:) = soilclass_default(:)
3300	> clayfraction(ib) = clayfraction_default
3301	> ELSE
3302	3571a3380,3381
3303	> DO ilf = 1,fopt
3304	> solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3305	3574c3384
3306	< ENDIF ! lucky/not lucky
3307	---
3308	> sgn = zero
3309	3576,3578c3386
3310	< ENDIF ! in the domain
3311	< ENDDO
3312	<
3313	---
3314	> ! Compute the average bare soil albedo parameters
3315	3580,3581c3388
3316	< ! Now we know how many points of which Olson type from the fine grid fall
3317	< ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3318	---
3319	> DO ilf = 1,fopt
3320	3583c3390
3321	<
3322	---
3323	> ! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean
3324	3585,3586c3392,3419
3325	< ! determine number of points of the fine grid which fall into each box of the
3326	< ! coarse grid
3327	---
3328	> IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND.&
3329	> & (solt(ilf) .NE. 6)) THEN
3330	> SELECTCASE(solt(ilf))
3331	> CASE(1)
3332	> soilclass(ib,1) = soilclass(ib,1) + sub_area(ib,ilf)
3333	> CASE(2)
3334	> soilclass(ib,2) = soilclass(ib,2) + sub_area(ib,ilf)
3335	> CASE(3)
3336	> soilclass(ib,2) = soilclass(ib,2) + sub_area(ib,ilf)
3337	> CASE(4)
3338	> soilclass(ib,2) = soilclass(ib,2) + sub_area(ib,ilf)
3339	> CASE(5)
3340	> soilclass(ib,3) = soilclass(ib,3) + sub_area(ib,ilf)
3341	> CASE(7)
3342	> soilclass(ib,2) = soilclass(ib,2) + sub_area(ib,ilf)
3343	> CASE DEFAULT
3344	> WRITE(numout,*) 'We should not be here, an impossible case appeared'
3345	> STOP 'slowproc_soilt'
3346	> END SELECT
3347	> clayfraction(ib) = clayfraction(ib) + &
3348	> & textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
3349	> sgn = sgn + sub_area(ib,ilf)
3350	> ELSE
3351	> IF (solt(ilf) .GT. nzobler) THEN
3352	> WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
3353	> STOP 'slowproc_soilt'
3354	> ENDIF
3355	> ENDIF
3356	3588,3589d3420
3357	< DO ib = 1, nbpt
3358	< n_found(ib) = SUM( n_origveg(ib,:) )
3359	3591,3593d3421
3360	<
3361	< !
3362	< ! determine fraction of Olson vegetation type in each box of the coarse grid
3363	3595,3602c3423
3364	< DO vid = 1, nolson
3365	< WHERE ( n_found(:) .GT. 0 )
3366	< frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std)
3367	< ELSEWHERE
3368	< frac_origveg(:,vid) = 0.
3369	< ENDWHERE
3370	< ENDDO
3371	<
3372	---
3373	> ! Normalize the surface
3374	3604,3605c3425,3432
3375	< ! now finally calculate coarse vegetation map
3376	< ! Find which model vegetation corresponds to each Olson type
3377	---
3378	> IF ( sgn .LT. min_sechiba) THEN
3379	> nbexp = nbexp + 1
3380	> soilclass(ib,:) = soilclass_default(:)
3381	> clayfraction(ib) = clayfraction_default
3382	> ELSE
3383	> soilclass(ib,:) = soilclass(ib,:)/sgn
3384	> clayfraction(ib) = clayfraction(ib)/sgn
3385	> ENDIF
3386	3607c3434
3387	< DO vid = 1, nolson
3388	---
3389	> ENDIF
3390	3609,3610d3435
3391	< DO jv = 1, nvm
3392	< veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3393	3611a3437
3394	>
3395	3613,3615c3439
3396	< DO jv = 1, nnobio
3397	< frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3398	< ENDDO
3399	---
3400	> CASE("fao")
3401	3617c3441
3402	< ENDDO
3403	---
3404	> soilclass_default=soilclass_default_fao
3405	3618a3443
3406	> PRINT *, "Using a soilclass map with fao classification"
3407	3620,3621c3445
3408	< WRITE(numout,*)
3409	< WRITE(numout,*) 'Interpolation Done'
3410	---
3411	> ALLOCATE(textfrac_table(nscm,ntext))
3412	3623c3447
3413	< ! Clean up the point of the map
3414	---
3415	> CALL get_soilcorr (nscm, textfrac_table)
3416	3627,3629c3451
3417	< ! Let us see if all points found something in the 5km map !
3418	< !
3419	< IF ( n_found(ib) .EQ. 0 ) THEN
3420	---
3421	> ! GO through the point we have found
3422	3631d3452
3423	< ! Now we need to handle some exceptions
3424	3633,3637c3454
3425	< IF ( lalo(ib,1) .LT. -56.0) THEN
3426	< ! Antartica
3427	< frac_nobio(ib,:) = 0.0
3428	< frac_nobio(ib,iice) = 1.0
3429	< veget(ib,:) = 0.0
3430	---
3431	> fopt = COUNT(sub_area(ib,:) > zero)
3432	3639,3643c3456
3433	< ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3434	< ! Artica
3435	< frac_nobio(ib,:) = 0.0
3436	< frac_nobio(ib,iice) = 1.0
3437	< veget(ib,:) = 0.0
3438	---
3439	> ! Check that we found some points
3440	3645,3649c3458,3459
3441	< ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3442	< ! Greenland
3443	< frac_nobio(ib,:) = 0.0
3444	< frac_nobio(ib,iice) = 1.0
3445	< veget(ib,:) = 0.0
3446	---
3447	> soilclass(ib,:) = 0.0
3448	> clayfraction(ib) = 0.0
3449	3650a3461,3464
3450	> IF ( fopt .EQ. 0) THEN
3451	> nbexp = nbexp + 1
3452	> soilclass(ib,:) = soilclass_default(:)
3453	> clayfraction(ib) = clayfraction_default
3454	3653,3655c3467,3469
3455	< WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box'
3456	< WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib)
3457	< WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib)
3458	---
3459	> DO ilf = 1,fopt
3460	> solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3461	> ENDDO
3462	3657,3661d3470
3463	< WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3464	< CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3465	< lalo(ib,2), lalo(ib,1), inear)
3466	< WRITE(numout,*) 'Coordinates of the nearest point:', &
3467	< lon_ful(inear),lat_ful(inear)
3468	3663,3665c3472
3469	< DO jv = 1, nvm
3470	< veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3471	< ENDDO
3472	---
3473	> ! Compute the average bare soil albedo parameters
3474	3667,3669c3474,3476
3475	< DO jv = 1, nnobio
3476	< frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3477	< ENDDO
3478	---
3479	> sgn = zero
3480	> !
3481	> DO ilf = 1,fopt
3482	3671d3477
3483	< ENDIF
3484	3672a3479,3488
3485	> !
3486	> IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3487	> soilclass(ib,solt(ilf)) = soilclass(ib,solt(ilf)) + sub_area(ib,ilf)
3488	> clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
3489	> sgn = sgn + sub_area(ib,ilf)
3490	> ELSE
3491	> IF (solt(ilf) .GT. nscm) THEN
3492	> WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3493	> STOP 'slowproc_soilt'
3494	> ENDIF
3495	3674a3491
3496	> ENDDO
3497	3676c3493
3498	< ! Limit the smalest vegetation fraction to 0.5%
3499	---
3500	> ! Normalize the surface
3501	3678,3680c3495,3501
3502	< DO vid = 1, nvm
3503	< IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3504	< veget(ib,vid) = 0.0
3505	---
3506	> IF ( sgn .LT. min_sechiba) THEN
3507	> nbexp = nbexp + 1
3508	> soilclass(ib,:) = soilclass_default(:)
3509	> clayfraction(ib) = clayfraction_default
3510	> ELSE
3511	> soilclass(ib,:) = soilclass(ib,:)/sgn
3512	> clayfraction(ib) = clayfraction(ib)/sgn
3513	3682,3686d3502
3514	< ENDDO
3515	< !
3516	< sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3517	< frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3518	< veget(ib,:) = veget(ib,:)/sumf
3519	3687a3504
3520	> ENDIF
3521	3691,3697d3507
3522	< DEALLOCATE(lon_up)
3523	< DEALLOCATE(lon_low)
3524	< DEALLOCATE(lat_up)
3525	< DEALLOCATE(lat_low)
3526	< DEALLOCATE(lat_ful)
3527	< DEALLOCATE(lon_ful)
3528	< DEALLOCATE(vegmap)
3529	3699c3509
3530	< RETURN
3531	---
3532	> CASE("fao2")
3533	3701,3705c3511
3534	< END SUBROUTINE slowproc_interpol_OLD_g
3535	< !!
3536	< !! Interpolate the IGBP vegetation map to the grid of the model
3537	< !!
3538	< SUBROUTINE slowproc_interpol_NEW_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
3539	---
3540	> soilclass_default=soilclass_default_fao2
3541	3706a3513
3542	> PRINT *, "Using a soilclass map with fao2 classification"
3543	3707a3515
3544	> ALLOCATE(textfrac_table(nscm,ntext))
3545	3709c3517
3546	< ! 0.1 INPUT
3547	---
3548	> CALL get_soilcorr (nscm, textfrac_table)
3549	3711,3727c3519
3550	< INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3551	< REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3552	< INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
3553	< ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
3554	< REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3555	< REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
3556	< !
3557	< ! 0.2 OUTPUT
3558	< !
3559	< REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
3560	< REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
3561	< !
3562	< LOGICAL :: ok_interpol ! optionnal return of aggregate_vec
3563	< !
3564	< ! 0.3 LOCAL
3565	< !
3566	< INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
3567	---
3568	> DO ib =1, nbpt
3569	3728a3521
3570	> ! GO through the point we have found
3571	3730,3742d3522
3572	< CHARACTER(LEN=80) :: filename
3573	< INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid
3574	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
3575	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3576	< INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index
3577	< REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg
3578	< REAL(r_std), DIMENSION(nbpt) :: n_found
3579	< REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
3580	< REAL(r_std) :: vegcorr(nolson,nvm)
3581	< REAL(r_std) :: nobiocorr(nolson,nnobio)
3582	< CHARACTER(LEN=40) :: callsign
3583	< REAL(r_std) :: sumf, resol_lon, resol_lat
3584	< INTEGER(i_std) :: idi, jv, inear, nbvmax
3585	3744c3524
3586	< INTEGER :: ALLOC_ERR
3587	---
3588	> fopt = COUNT(sub_area(ib,:) > zero)
3589	3746,3747c3526
3590	< n_origveg(:,:) = zero
3591	< n_found(:) = zero
3592	---
3593	> ! Check that we found some points
3594	3749c3528,3529
3595	< CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3596	---
3597	> soilclass(ib,:) = 0.0
3598	> clayfraction(ib) = 0.0
3599	3751,3759c3531,3535
3600	< !Config Key = VEGETATION_FILE
3601	< !Config Desc = Name of file from which the vegetation map is to be read
3602	< !Config If = !IMPOSE_VEG
3603	< !Config If = !LAND_USE
3604	< !Config Def = ../surfmap/carteveg5km.nc
3605	< !Config Help = The name of the file to be opened to read the vegetation
3606	< !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3607	< !Config map which is derived from the IGBP one. We assume that we have
3608	< !Config a classification in 87 types. This is Olson modified by Viovy.
3609	---
3610	> IF ( fopt .EQ. 0) THEN
3611	> nbexp = nbexp + 1
3612	> soilclass(ib,:) = soilclass_default(:)
3613	> clayfraction(ib) = clayfraction_default
3614	> ELSE
3615	3761,3762c3537,3540
3616	< filename = '../surfmap/carteveg5km.nc'
3617	< CALL getin('VEGETATION_FILE',filename)
3618	---
3619	> DO ilf = 1,fopt
3620	> solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3621	> solt2(ilf) = soiltext2(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3622	> ENDDO
3623	3764d3541
3624	< CALL flininfo(filename, iml, jml, lml, tml, fid)
3625	3765a3543
3626	> ! Compute the average bare soil albedo parameters
3627	3767,3784c3545
3628	< ALLOC_ERR=-1
3629	< ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR)
3630	< IF (ALLOC_ERR/=0) THEN
3631	< WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
3632	< STOP
3633	< ENDIF
3634	< ALLOC_ERR=-1
3635	< ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR)
3636	< IF (ALLOC_ERR/=0) THEN
3637	< WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
3638	< STOP
3639	< ENDIF
3640	< ALLOC_ERR=-1
3641	< ALLOCATE(vegmap(iml), STAT=ALLOC_ERR)
3642	< IF (ALLOC_ERR/=0) THEN
3643	< WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
3644	< STOP
3645	< ENDIF
3646	---
3647	> sgn = zero
3648	3786c3547
3649	< WRITE(numout,*) 'Reading the OLSON type vegetation file'
3650	---
3651	> DO ilf = 1,fopt
3652	3788,3790d3548
3653	< CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
3654	< CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
3655	< CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
3656	3792d3549
3657	< CALL flinclo(fid)
3658	3794,3801c3551,3555
3659	< IF (MAXVAL(vegmap) .LT. nolson) THEN
3660	< WRITE(numout,*) 'WARNING -- WARNING'
3661	< WRITE(numout,*) 'The vegetation map has to few vegetation types.'
3662	< WRITE(numout,*) 'If you are lucky it will work but please check'
3663	< ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
3664	< WRITE(numout,*) 'More vegetation types in file than the code can'
3665	< WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson
3666	< STOP 'slowproc_interpol'
3667	---
3668	> IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3669	> IF ( solt2(ilf) .GT. solt(ilf)) THEN
3670	> soilclass(ib,2solt(ilf)) = soilclass(ib,2solt(ilf)) + sub_area(ib,ilf)
3671	> ELSE
3672	> soilclass(ib,2solt(ilf)-1) = soilclass(ib,2solt(ilf)-1) + sub_area(ib,ilf)
3673	3803,3826c3557,3562
3674	< !
3675	< ! Some assumptions on the vegetation file. This information should be
3676	< ! be computed or read from the file.
3677	< ! It is the reolution in meters of the grid of the vegetation file.
3678	< !
3679	< resol_lon = 5000.
3680	< resol_lat = 5000.
3681	< !
3682	< ! The number of maximum vegetation map points in the GCM grid should
3683	< ! also be computed and not imposed here.
3684	< nbvmax = iml/nbpt
3685	< !
3686	< callsign="Vegetation map"
3687	< !
3688	< ok_interpol = .FALSE.
3689	< DO WHILE ( .NOT. ok_interpol )
3690	< WRITE(numout,*) "Projection arrays for ",callsign," : "
3691	< WRITE(numout,*) "nbvmax = ",nbvmax
3692	< !
3693	< ALLOC_ERR=-1
3694	< ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
3695	< IF (ALLOC_ERR/=0) THEN
3696	< WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3697	< STOP
3698	---
3699	> clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
3700	> sgn = sgn + sub_area(ib,ilf)
3701	> ELSE
3702	> IF (solt(ilf) .GT. nscm) THEN
3703	> WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3704	> STOP 'slowproc_soilt'
3705	3828,3833d3563
3706	< sub_index(:,:)=0
3707	< ALLOC_ERR=-1
3708	< ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
3709	< IF (ALLOC_ERR/=0) THEN
3710	< WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3711	< STOP
3712	3835d3564
3713	< sub_area(:,:)=zero
3714	3837,3839c3566
3715	< CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, &
3716	< & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
3717	< & nbvmax, sub_index, sub_area, ok_interpol)
3718	---
3719	> ENDDO
3720	3841,3843c3568
3721	< IF ( .NOT. ok_interpol ) THEN
3722	< DEALLOCATE(sub_area)
3723	< DEALLOCATE(sub_index)
3724	---
3725	> ! Normalize the surface
3726	3845c3570,3573
3727	< nbvmax = nbvmax * 2
3728	---
3729	> IF ( sgn .LT. min_sechiba) THEN
3730	> nbexp = nbexp + 1
3731	> soilclass(ib,:) = soilclass_default(:)
3732	> clayfraction(ib) = clayfraction_default
3733	3847,3857c3575,3576
3734	< !
3735	< DO ib = 1, nbpt
3736	< idi=1
3737	< DO WHILE ( sub_area(ib,idi) > zero )
3738	< ip = sub_index(ib,idi)
3739	< n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi)
3740	< n_found(ib) = n_found(ib) + sub_area(ib,idi)
3741	< idi = idi +1
3742	< ENDDO
3743	< ENDDO
3744	< !
3745	---
3746	> soilclass(ib,:) = soilclass(ib,:)/sgn
3747	> clayfraction(ib) = clayfraction(ib)/sgn
3748	3859,3863d3577
3749	< ENDDO
3750	< !
3751	< ! Now we know how many points of which Olson type from the fine grid fall
3752	< ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3753	< !
3754	3865c3579
3755	< ! determine fraction of Olson vegetation type in each box of the coarse grid
3756	---
3757	> ENDIF
3758	3867,3872d3580
3759	< DO vid = 1, nolson
3760	< WHERE ( n_found(:) .GT. 0 )
3761	< frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:)
3762	< ELSEWHERE
3763	< frac_origveg(:,vid) = 0.
3764	< ENDWHERE
3765	3875,3876d3582
3766	< ! now finally calculate coarse vegetation map
3767	< ! Find which model vegetation corresponds to each Olson type
3768	3878,3879c3584
3769	< veget(:,:) = zero
3770	< frac_nobio(:,:) = zero
3771	---
3772	> CASE("usda")
3773	3881c3586
3774	< DO vid = 1, nolson
3775	---
3776	> soilclass_default=soilclass_default_usda
3777	3883,3891c3588
3778	< DO jv = 1, nvm
3779	< veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3780	< ENDDO
3781	< !
3782	< DO jv = 1, nnobio
3783	< frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3784	< ENDDO
3785	< !
3786	< ENDDO
3787	---
3788	> PRINT *, "Using a soilclass map with usda classification"
3789	3893c3590
3790	< WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
3791	---
3792	> ALLOCATE(textfrac_table(nscm,ntext))
3793	3895c3592
3794	< ! Clean up the point of the map
3795	---
3796	> CALL get_soilcorr (nscm, textfrac_table)
3797	3899c3596
3798	< ! Let us see if all points found something in the 5km map !
3799	---
3800	> ! GO through the point we have found
3801	3901d3597
3802	< IF ( n_found(ib) .EQ. 0 ) THEN
3803	3903c3599
3804	< ! Now we need to handle some exceptions
3805	---
3806	> fopt = COUNT(sub_area(ib,:) > zero)
3807	3905,3909c3601
3808	< IF ( lalo(ib,1) .LT. -56.0) THEN
3809	< ! Antartica
3810	< frac_nobio(ib,:) = 0.0
3811	< frac_nobio(ib,iice) = 1.0
3812	< veget(ib,:) = 0.0
3813	---
3814	> ! Check that we found some points
3815	3911,3921c3603,3604
3816	< ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3817	< ! Artica
3818	< frac_nobio(ib,:) = 0.0
3819	< frac_nobio(ib,iice) = 1.0
3820	< veget(ib,:) = 0.0
3821	< !
3822	< ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3823	< ! Greenland
3824	< frac_nobio(ib,:) = 0.0
3825	< frac_nobio(ib,iice) = 1.0
3826	< veget(ib,:) = 0.0
3827	---
3828	> soilclass(ib,:) = 0.0
3829	> clayfraction(ib) = 0.0
3830	3922a3606,3609
3831	> IF ( fopt .EQ. 0) THEN
3832	> nbexp = nbexp + 1
3833	> soilclass(ib,:) = soilclass_default
3834	> clayfraction(ib) = clayfraction_default
3835	3925,3940c3612,3613
3836	< WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib
3837	< WRITE(numout,*) 'Longitude range : ', lalo(ib,2)
3838	< WRITE(numout,*) 'Latitude range : ', lalo(ib,1)
3839	< !
3840	< WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3841	< CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3842	< lalo(ib,2), lalo(ib,1), inear)
3843	< WRITE(numout,*) 'Coordinates of the nearest point:', &
3844	< lon_ful(inear),lat_ful(inear)
3845	< !
3846	< DO jv = 1, nvm
3847	< veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3848	< ENDDO
3849	< !
3850	< DO jv = 1, nnobio
3851	< frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3852	---
3853	> DO ilf = 1,fopt
3854	> solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3855	3943d3615
3856	< ENDIF
3857	3945c3617
3858	< ENDIF
3859	---
3860	> ! Compute the average bare soil albedo parameters
3861	3946a3619
3862	> sgn = zero
3863	3948c3621
3864	< ! Limit the smalest vegetation fraction to 0.5%
3865	---
3866	> DO ilf = 1,fopt
3867	3950,3954d3622
3868	< DO vid = 1, nvm
3869	< IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3870	< veget(ib,vid) = 0.0
3871	< ENDIF
3872	< ENDDO
3873	3956,3958d3623
3874	< sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3875	< frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3876	< veget(ib,:) = veget(ib,:)/sumf
3877	3959a3625,3634
3878	> IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3879	> soilclass(ib,solt(ilf)) = soilclass(ib,solt(ilf)) + sub_area(ib,ilf)
3880	> clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
3881	> sgn = sgn + sub_area(ib,ilf)
3882	> ELSE
3883	> IF (solt(ilf) .GT. nscm) THEN
3884	> WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3885	> STOP 'slowproc_soilt'
3886	> ENDIF
3887	> ENDIF
3888	3963,4016c3638,3650
3889	< DEALLOCATE(vegmap)
3890	< DEALLOCATE(lat_ful, lon_ful)
3891	< DEALLOCATE(sub_index)
3892	< DEALLOCATE(sub_area)
3893	<
3894	< !
3895	< RETURN
3896	< !
3897	< END SUBROUTINE slowproc_interpol_NEW_g
3898	<
3899	<
3900	< !!
3901	< !! looks for nearest grid point on the fine map
3902	< !!
3903	< SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
3904	<
3905	< INTEGER(i_std), INTENT(in) :: iml
3906	< REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5
3907	< REAL(r_std), INTENT(in) :: lonmod, latmod
3908	<
3909	< INTEGER(i_std), INTENT(out) :: inear
3910	<
3911	< REAL(r_std) :: pi
3912	< REAL(r_std) :: pa, p
3913	< REAL(r_std) :: coscolat, sincolat
3914	< REAL(r_std) :: cospa, sinpa
3915	< REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
3916	< INTEGER(i_std) :: i
3917	< INTEGER(i_std), DIMENSION(1) :: ineartab
3918	< INTEGER :: ALLOC_ERR
3919	<
3920	< ALLOC_ERR=-1
3921	< ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
3922	< IF (ALLOC_ERR/=0) THEN
3923	< WRITE(numout,*) "ERROR IN ALLOCATION of cosang : ",ALLOC_ERR
3924	< STOP
3925	< ENDIF
3926	<
3927	< pi = 4.0 * ATAN(1.0)
3928	<
3929	< pa = pi/2.0 - latmod*pi/180.0 ! dist. entre pole n et point a
3930	< cospa = COS(pa)
3931	< sinpa = SIN(pa)
3932	<
3933	< DO i = 1, iml
3934	<
3935	< sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 )
3936	< coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 )
3937	<
3938	< p = (lonmod-lon5(i))*pi/180.0 ! angle entre a et b (leurs meridiens)
3939	<
3940	< ! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
3941	< cosang(i) = cospacoscolat + sinpasincolat*COS(p)
3942	<
3943	---
3944	> ! Normalize the surface
3945	> !
3946	> IF ( sgn .LT. min_sechiba) THEN
3947	> nbexp = nbexp + 1
3948	> soilclass(ib,:) = soilclass_default(:)
3949	> clayfraction(ib) = clayfraction_default
3950	> ELSE
3951	> soilclass(ib,:) = soilclass(ib,:)/sgn
3952	> clayfraction(ib) = clayfraction(ib)/sgn
3953	> ENDIF
3954	> !
3955	> ENDIF
3956	> !
3957	4018,4023c3652,3685
3958	<
3959	< ineartab = MAXLOC( cosang(:) )
3960	< inear = ineartab(1)
3961	<
3962	< DEALLOCATE(cosang)
3963	< END SUBROUTINE slowproc_nearest
3964	---
3965	> !
3966	> CASE DEFAULT
3967	> WRITE(,) 'A non supported soil type classification has been chosen'
3968	> STOP 'slowproc_soilt'
3969	> ENDSELECT
3970	> !
3971	> WRITE(numout,*) 'Interpolation Done'
3972	> !
3973	> IF ( nbexp .GT. 0 ) THEN
3974	> WRITE(numout,*) 'slowproc_soilt : The interpolation of the bare soil albedo had ', nbexp
3975	> WRITE(numout,*) 'slowproc_soilt : points without data. This are either coastal points or'
3976	> WRITE(numout,*) 'slowproc_soilt : ice covered land.'
3977	> WRITE(numout,*) 'slowproc_soilt : The problem was solved by using the default soil types.'
3978	> ENDIF
3979	> !
3980	> DEALLOCATE (lat_rel)
3981	> DEALLOCATE (lon_rel)
3982	> DEALLOCATE (mask)
3983	> DEALLOCATE (sub_area)
3984	> DEALLOCATE (sub_index)
3985	> DEALLOCATE (soiltext)
3986	> DEALLOCATE (solt)
3987	> !
3988	> DEALLOCATE (soiltext2)
3989	> DEALLOCATE (solt2)
3990	> DEALLOCATE (textfrac_table)
3991	> !
3992	> DEALLOCATE (resol_lu)
3993	> !
3994	> !
3995	> RETURN
3996	> !
3997	> END SUBROUTINE slowproc_soilt
3998	> !
3999	4026c3688
4000	< !! Interpolate the Zobler soil type map
4001	---
4002	> !! Calculate mean slope coef in each model grid box from the slope map
4003	4028,4029c3690
4004	< SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soiltype, clayfraction)
4005	< !
4006	---
4007	> SUBROUTINE slowproc_slope(nbpt, lalo, neighbours, resolution, contfrac, reinf_slope)
4008	4031,4034d3691
4009	< ! This subroutine should read the Zobler map and interpolate to the model grid. The method
4010	< ! is to get fraction of the three main soiltypes for each grid box.
4011	< ! The soil fraction are going to be put into the array soiltype in the following order :
4012	< ! coarse, medium and fine.
4013	4044c3701
4014	< REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
4015	---
4016	> REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Fraction of continent in the grid
4017	4048,4050c3705
4018	< REAL(r_std), INTENT(out) :: soiltype(nbpt, nstm) ! Soil type map to be created from the Zobler map
4019	< REAL(r_std), INTENT(out) :: clayfraction(nbpt) ! The fraction of clay as used by STOMATE
4020	< !
4021	---
4022	> REAL(r_std), INTENT(out) :: reinf_slope(nbpt) ! slope coef
4023	4054d3708
4024	< INTEGER(i_std) :: nbvmax
4025	4055a3710
4026	> REAL(r_std) :: slope_noreinf ! Slope above which runoff is maximum
4027	4057,4060c3712,3715
4028	< INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, nbexp
4029	< REAL(r_std) :: lev(1), date, dt
4030	< INTEGER(i_std) :: itau(1)
4031	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, soiltext
4032	---
4033	> CHARACTER(LEN=30) :: callsign
4034	> INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid
4035	> INTEGER(i_std) :: idi, nbvmax
4036	> REAL(r_std) :: slopecoef, coslat
4037	4062d3716
4038	< REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
4039	4064,4071c3718,3722
4040	< INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt
4041	< REAL(r_std) :: sgn
4042	< CHARACTER(LEN=30) :: callsign
4043	< !
4044	< ! Number of texture classes in Zobler
4045	< !
4046	< INTEGER(i_std), PARAMETER :: classnb = 7
4047	< REAL(r_std) :: textfrac_table(classnb, nstm)
4048	---
4049	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, slopemap
4050	> REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu
4051	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
4052	> REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu
4053	> INTEGER(i_std) :: nix, njx
4054	4073c3724
4055	< LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
4056	---
4057	> LOGICAL :: ok_interpol = .FALSE. ! optionnal return of aggregate_2d
4058	4078d3728
4059	< CALL get_soilcorr (classnb, textfrac_table)
4060	4080c3730,3734
4061	< ! Needs to be a configurable variable
4062	---
4063	> !Config Key = SLOPE_NOREINF
4064	> !Config Desc = See slope_noreinf above
4065	> !Config If =
4066	> !Config Def = 0.5
4067	> !Config Help = The slope above which there is no reinfiltration
4068	4081a3736
4069	> slope_noreinf = 0.5
4070	4083,4091c3738
4071	< !Config Key = SOILTYPE_FILE
4072	< !Config Desc = Name of file from which soil types are read
4073	< !Config Def = ../surfmap/soils_param.nc
4074	< !Config If = !IMPOSE_VEG
4075	< !Config Help = The name of the file to be opened to read the soil types.
4076	< !Config The data from this file is then interpolated to the grid of
4077	< !Config of the model. The aim is to get fractions for sand loam and
4078	< !Config clay in each grid box. This information is used for soil hydrology
4079	< !Config and respiration.
4080	---
4081	> CALL getin_p('SLOPE_NOREINF',slope_noreinf)
4082	4093,4094c3740,3746
4083	< filename = '../surfmap/soils_param.nc'
4084	< CALL getin_p('SOILTYPE_FILE',filename)
4085	---
4086	> !Config Key = TOPOGRAPHY_FILE
4087	> !Config Desc = Name of file from which the topography map is to be read
4088	> !Config If =
4089	> !Config Def = ../surfmap/cartepente.nc
4090	> !Config Help = The name of the file to be opened to read the orography
4091	> !Config map is to be given here. Usualy SECHIBA runs with a 2'
4092	> !Config map which is derived from the NGDC one.
4093	4096,4099c3748,3751
4094	< IF (is_root_prc) THEN
4095	< CALL flininfo(filename,iml, jml, lml, tml, fid)
4096	< CALL flinclo(fid)
4097	< ENDIF
4098	---
4099	> filename = '../surfmap/cartepente2d_15min.nc'
4100	> CALL getin_p('TOPOGRAPHY_FILE',filename)
4101	> !
4102	> IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
4103	4105,4106d3756
4104	< ! soils_param.nc file is 1° soit texture file.
4105	< !
4106	4108c3758
4107	< ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
4108	---
4109	> ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
4110	4110c3760
4111	< WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
4112	---
4113	> WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
4114	4114c3764
4115	< ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
4116	---
4117	> ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
4118	4116c3766
4119	< WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
4120	---
4121	> WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
4122	4120c3770
4123	< ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4124	---
4125	> ALLOCATE(slopemap(iml,jml), STAT=ALLOC_ERR)
4126	4122c3772
4127	< WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4128	---
4129	> WRITE(numout,*) "ERROR IN ALLOCATION of slopemap : ",ALLOC_ERR
4130	4126c3776
4131	< ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR)
4132	---
4133	> ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
4134	4128c3778
4135	< WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
4136	---
4137	> WRITE(numout,*) "ERROR IN ALLOCATION of resol_lu : ",ALLOC_ERR
4138	4132,4138c3782
4139	< IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
4140	< CALL bcast(lon_rel)
4141	< CALL bcast(lat_rel)
4142	< CALL bcast(itau)
4143	< CALL bcast(date)
4144	< CALL bcast(dt)
4145	<
4146	---
4147	> WRITE(numout,*) 'Reading the topography file'
4148	4140,4141c3784,3787
4149	< IF (is_root_prc) CALL flinget(fid, 'soiltext', iml, jml, lml, tml, 1, 1, soiltext)
4150	< CALL bcast(soiltext)
4151	---
4152	> IF (is_root_prc) THEN
4153	> CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
4154	> CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
4155	> CALL flinget(fid, 'pente', iml, jml, 0, 0, 1, 1, slopemap)
4156	4143c3789,3793
4157	< IF (is_root_prc) CALL flinclo(fid)
4158	---
4159	> CALL flinclo(fid)
4160	> ENDIF
4161	> CALL bcast(lon_lu)
4162	> CALL bcast(lat_lu)
4163	> CALL bcast(slopemap)
4164	4145c3795,3813
4165	< nbexp = 0
4166	---
4167	> ALLOC_ERR=-1
4168	> ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
4169	> IF (ALLOC_ERR/=0) THEN
4170	> WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
4171	> STOP
4172	> ENDIF
4173	> ALLOC_ERR=-1
4174	> ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
4175	> IF (ALLOC_ERR/=0) THEN
4176	> WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
4177	> STOP
4178	> ENDIF
4179	> !
4180	> DO ip=1,iml
4181	> lat_rel(ip,:) = lat_lu(:)
4182	> ENDDO
4183	> DO jp=1,jml
4184	> lon_rel(:,jp) = lon_lu(:)
4185	> ENDDO
4186	4149a3818,3824
4187	> ALLOC_ERR=-1
4188	> ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4189	> IF (ALLOC_ERR/=0) THEN
4190	> WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4191	> STOP
4192	> ENDIF
4193	> !
4194	4153c3828
4195	< IF (soiltext(ip,jp) .GT. min_sechiba) THEN
4196	---
4197	> IF (slopemap(ip,jp) .GT. min_sechiba) THEN
4198	4155a3831,3852
4199	> !
4200	> ! Resolution in longitude
4201	> !
4202	> coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )
4203	> IF ( ip .EQ. 1 ) THEN
4204	> resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip,jp) ) * pi/180. * R_Earth * coslat
4205	> ELSEIF ( ip .EQ. iml ) THEN
4206	> resol_lu(ip,jp,1) = ABS( lon_rel(ip,jp) - lon_rel(ip-1,jp) ) * pi/180. * R_Earth * coslat
4207	> ELSE
4208	> resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip-1,jp) )/2. * pi/180. * R_Earth * coslat
4209	> ENDIF
4210	> !
4211	> ! Resolution in latitude
4212	> !
4213	> IF ( jp .EQ. 1 ) THEN
4214	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp) - lat_rel(ip,jp+1) ) * pi/180. * R_Earth
4215	> ELSEIF ( jp .EQ. jml ) THEN
4216	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp-1) - lat_rel(ip,jp) ) * pi/180. * R_Earth
4217	> ELSE
4218	> resol_lu(ip,jp,2) = ABS( lat_rel(ip,jp-1) - lat_rel(ip,jp+1) )/2. * pi/180. * R_Earth
4219	> ENDIF
4220	> !
4221	4159d3855
4222	< nbvmax = 200
4223	4161c3857,3867
4224	< callsign = "Soil types"
4225	---
4226	> ! The number of maximum vegetation map points in the GCM grid is estimated.
4227	> ! Some lmargin is taken.
4228	> !
4229	> IF (is_root_prc) THEN
4230	> nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2
4231	> njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2
4232	> nbvmax = nix*njx
4233	> ENDIF
4234	> CALL bcast(nbvmax)
4235	> !
4236	> callsign="Slope map"
4237	4163,4164d3868
4238	< ok_interpol = .FALSE.
4239	< DO WHILE ( .NOT. ok_interpol )
4240	4169,4174d3872
4241	< ALLOCATE(solt(nbvmax), STAT=ALLOC_ERR)
4242	< IF (ALLOC_ERR/=0) THEN
4243	< WRITE(numout,*) "ERROR IN ALLOCATION of solt : ",ALLOC_ERR
4244	< STOP
4245	< ENDIF
4246	< ALLOC_ERR=-1
4247	4193,4200d3890
4248	< IF ( .NOT. ok_interpol ) THEN
4249	< DEALLOCATE(sub_area)
4250	< DEALLOCATE(sub_index)
4251	< DEALLOCATE(solt)
4252	< !
4253	< nbvmax = nbvmax * 2
4254	< ENDIF
4255	< ENDDO
4256	4202a3893,3900
4257	> idi=1
4258	> !-
4259	> !- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro
4260	> !-
4261	> slopecoef = zero
4262	> DO WHILE ( sub_area(ib,idi) > zero )
4263	> ip = sub_index(ib,idi,1)
4264	> jp = sub_index(ib,idi,2)
4265	4204,4221c3902,3903
4266	< soiltype(ib,:) = zero
4267	< clayfraction(ib) = zero
4268	< !
4269	< ! GO through the point we have found
4270	< !
4271	< !
4272	< fopt = COUNT(sub_area(ib,:) > zero)
4273	< !
4274	< ! Check that we found some points
4275	< !
4276	< IF ( fopt .EQ. 0) THEN
4277	< nbexp = nbexp + 1
4278	< soiltype(ib,:) = soiltype_default(:)
4279	< clayfraction(ib) = clayfraction_default
4280	< ELSE
4281	< !
4282	< DO ilf = 1,fopt
4283	< solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
4284	---
4285	> slopecoef = slopecoef + MIN(slopemap(ip,jp)/slope_noreinf, un) * sub_area(ib,idi)
4286	> idi = idi +1
4287	4223,4253c3905,3906
4288	< !
4289	< sgn = zero
4290	< !
4291	< ! Compute the average bare soil albedo parameters
4292	< !
4293	< DO ilf = 1,fopt
4294	< !
4295	< ! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean
4296	< !
4297	< IF ( (solt(ilf) .LE. classnb) .AND. (solt(ilf) .GT. 0) .AND.&
4298	< & (solt(ilf) .NE. 6)) THEN
4299	< SELECTCASE(solt(ilf))
4300	< CASE(1)
4301	< soiltype(ib,1) = soiltype(ib,1) + sub_area(ib,ilf)
4302	< CASE(2)
4303	< soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4304	< CASE(3)
4305	< soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4306	< CASE(4)
4307	< soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4308	< CASE(5)
4309	< soiltype(ib,3) = soiltype(ib,3) + sub_area(ib,ilf)
4310	< CASE(7)
4311	< soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4312	< CASE DEFAULT
4313	< WRITE(numout,*) 'We should not be here, an impossible case appeared'
4314	< STOP 'slowproc_soilt'
4315	< END SELECT
4316	< clayfraction(ib) = clayfraction(ib) + &
4317	< & textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
4318	< sgn = sgn + sub_area(ib,ilf)
4319	---
4320	> IF ( idi .GT. 1 ) THEN
4321	> reinf_slope(ib) = un - slopecoef / SUM(sub_area(ib,1:(idi-1)))
4322	4255,4258c3908
4323	< IF (solt(ilf) .GT. classnb) THEN
4324	< WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
4325	< STOP 'slowproc_soilt'
4326	< ENDIF
4327	---
4328	> reinf_slope(ib) = slope_default
4329	4260d3909
4330	< !
4331	4263,4274d3911
4332	< ! Normalize the surface
4333	< !
4334	< IF ( sgn .LT. min_sechiba) THEN
4335	< nbexp = nbexp + 1
4336	< soiltype(ib,:) = soiltype_default(:)
4337	< clayfraction(ib) = clayfraction_default
4338	< ELSE
4339	< soiltype(ib,:) = soiltype(ib,:)/sgn
4340	< clayfraction(ib) = clayfraction(ib)/sgn
4341	< ENDIF
4342	< !
4343	< ENDIF
4344	4276c3913
4345	< ENDDO
4346	---
4347	> WRITE(numout,*) 'Interpolation Done'
4348	4278,4283d3914
4349	< IF ( nbexp .GT. 0 ) THEN
4350	< WRITE(numout,*) 'slowproc_soilt : The interpolation of the bare soil albedo had ', nbexp
4351	< WRITE(numout,*) 'slowproc_soilt : points without data. This are either coastal points or'
4352	< WRITE(numout,*) 'slowproc_soilt : ice covered land.'
4353	< WRITE(numout,*) 'slowproc_soilt : The problem was solved by using the default soil types.'
4354	< ENDIF
4355	4285,4288c3916
4356	< DEALLOCATE (lat_rel)
4357	< DEALLOCATE (lon_rel)
4358	< DEALLOCATE (mask)
4359	< DEALLOCATE (sub_area)
4360	---
4361	> DEALLOCATE(slopemap)
4362	4290,4291c3918,3923
4363	< DEALLOCATE (soiltext)
4364	< DEALLOCATE (solt)
4365	---
4366	> DEALLOCATE(sub_area)
4367	> DEALLOCATE(mask)
4368	> DEALLOCATE(lon_lu)
4369	> DEALLOCATE(lat_lu)
4370	> DEALLOCATE(lon_rel)
4371	> DEALLOCATE(lat_rel)
4372	4296,4298c3928
4373	< END SUBROUTINE slowproc_soilt
4374	< !
4375	<
4376	---
4377	> END SUBROUTINE slowproc_slope
4378	_______________________________________________________________________________________________________________________
4379	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE/src_sechiba/thermosoil.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba/thermosoil.f90
4380	107,109d106
4381	< REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp
4382	< REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1
4383	< REAL(r_std),DIMENSION (kjpindex) :: temp2
4384	129a127
4385	>
4386	135,229d132
4387	< IF (ldrestart_read) THEN
4388	< IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
4389	<
4390	< var_name= 'cgrnd'
4391	< CALL ioconf_setatt('UNITS', '-')
4392	< CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.')
4393	< IF ( ok_var(var_name) ) THEN
4394	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4395	< IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4396	< cgrnd(:,:)=temp1(:,:)
4397	< ENDIF
4398	< ENDIF
4399	<
4400	< var_name= 'dgrnd'
4401	< CALL ioconf_setatt('UNITS', '-')
4402	< CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.')
4403	< IF ( ok_var(var_name) ) THEN
4404	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4405	< IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4406	< dgrnd(:,:)=temp1(:,:)
4407	< ENDIF
4408	< ENDIF
4409	<
4410	< var_name= 'z1'
4411	< CALL ioconf_setatt('UNITS', '-')
4412	< CALL ioconf_setatt('LONG_NAME','?.')
4413	< IF ( ok_var(var_name) ) THEN
4414	< CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
4415	< IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
4416	< z1(:)=temp2(:)
4417	< ENDIF
4418	< ENDIF
4419	<
4420	< var_name= 'pcapa'
4421	< CALL ioconf_setatt('UNITS', '-')
4422	< CALL ioconf_setatt('LONG_NAME','?.')
4423	< IF ( ok_var(var_name) ) THEN
4424	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4425	< IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4426	< pcapa(:,:)=temp(:,:)
4427	< ENDIF
4428	< ENDIF
4429	<
4430	< var_name= 'pcapa_en'
4431	< CALL ioconf_setatt('UNITS', '-')
4432	< CALL ioconf_setatt('LONG_NAME','?.')
4433	< IF ( ok_var(var_name) ) THEN
4434	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4435	< IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4436	< pcapa_en(:,:)=temp(:,:)
4437	< ENDIF
4438	< ENDIF
4439	<
4440	< var_name= 'pkappa'
4441	< CALL ioconf_setatt('UNITS', '-')
4442	< CALL ioconf_setatt('LONG_NAME','?.')
4443	< IF ( ok_var(var_name) ) THEN
4444	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4445	< IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4446	< pkappa(:,:)=temp(:,:)
4447	< ENDIF
4448	< ENDIF
4449	<
4450	< var_name= 'zdz1'
4451	< CALL ioconf_setatt('UNITS', '-')
4452	< CALL ioconf_setatt('LONG_NAME','?.')
4453	< IF ( ok_var(var_name) ) THEN
4454	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4455	< IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4456	< zdz1(:,:)=temp1(:,:)
4457	< ENDIF
4458	< ENDIF
4459	<
4460	< var_name= 'zdz2'
4461	< CALL ioconf_setatt('UNITS', '-')
4462	< CALL ioconf_setatt('LONG_NAME','?.')
4463	< IF ( ok_var(var_name) ) THEN
4464	< CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4465	< IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4466	< zdz2(:,:)=temp(:,:)
4467	< ENDIF
4468	< ENDIF
4469	<
4470	< var_name='temp_sol_beg'
4471	< CALL ioconf_setatt('UNITS', 'K')
4472	< CALL ioconf_setatt('LONG_NAME','Old Surface temperature')
4473	< IF ( ok_var(var_name) ) THEN
4474	< CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
4475	< IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
4476	< temp_sol_beg(:) = temp2(:)
4477	< ENDIF
4478	< ENDIF
4479	<
4480	< ENDIF
4481	<
4482	318a222,224
4483	> !!$ WRITE(,) 'TTTTTX 2470 T :', mpi_rank, '--', temp_sol_new(2470), ptn(2470,1:3)
4484	> !!$ WRITE(,) 'TTTTTX 2470 F :', mpi_rank, '--', soilflx(2470), soilcap(2470)
4485	>
4486	341c247,251
4487	< INTEGER(i_std) :: ier
4488	---
4489	> INTEGER(i_std) :: ier, i_notfound
4490	> REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp
4491	> REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1
4492	> REAL(r_std),DIMENSION (kjpindex) :: temp2
4493	> REAL(r_std) :: ptn_default
4494	456c366
4495	<
4496	---
4497	> !
4498	463a374,485
4499	>
4500	> i_notfound = 0
4501	>
4502	> var_name= 'cgrnd'
4503	> CALL ioconf_setatt('UNITS', '-')
4504	> CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.')
4505	> IF ( ok_var(var_name) ) THEN
4506	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4507	> IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4508	> cgrnd(:,:)=temp1(:,:)
4509	> ELSE
4510	> i_notfound = i_notfound + 1
4511	> ENDIF
4512	> ENDIF
4513	>
4514	> var_name= 'dgrnd'
4515	> CALL ioconf_setatt('UNITS', '-')
4516	> CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.')
4517	> IF ( ok_var(var_name) ) THEN
4518	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4519	> IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4520	> dgrnd(:,:)=temp1(:,:)
4521	> ELSE
4522	> i_notfound = i_notfound + 1
4523	> ENDIF
4524	> ENDIF
4525	>
4526	> var_name= 'z1'
4527	> CALL ioconf_setatt('UNITS', '-')
4528	> CALL ioconf_setatt('LONG_NAME','?.')
4529	> IF ( ok_var(var_name) ) THEN
4530	> CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
4531	> IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
4532	> z1(:)=temp2(:)
4533	> ELSE
4534	> i_notfound = i_notfound + 1
4535	> ENDIF
4536	> ENDIF
4537	>
4538	> var_name= 'pcapa'
4539	> CALL ioconf_setatt('UNITS', '-')
4540	> CALL ioconf_setatt('LONG_NAME','?.')
4541	> IF ( ok_var(var_name) ) THEN
4542	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4543	> IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4544	> pcapa(:,:)=temp(:,:)
4545	> ELSE
4546	> i_notfound = i_notfound + 1
4547	> ENDIF
4548	> ENDIF
4549	>
4550	> var_name= 'pcapa_en'
4551	> CALL ioconf_setatt('UNITS', '-')
4552	> CALL ioconf_setatt('LONG_NAME','?.')
4553	> IF ( ok_var(var_name) ) THEN
4554	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4555	> IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4556	> pcapa_en(:,:)=temp(:,:)
4557	> ELSE
4558	> i_notfound = i_notfound + 1
4559	> ENDIF
4560	> ENDIF
4561	>
4562	> var_name= 'pkappa'
4563	> CALL ioconf_setatt('UNITS', '-')
4564	> CALL ioconf_setatt('LONG_NAME','?.')
4565	> IF ( ok_var(var_name) ) THEN
4566	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4567	> IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4568	> pkappa(:,:)=temp(:,:)
4569	> ELSE
4570	> i_notfound = i_notfound + 1
4571	> ENDIF
4572	> ENDIF
4573	>
4574	> var_name= 'zdz1'
4575	> CALL ioconf_setatt('UNITS', '-')
4576	> CALL ioconf_setatt('LONG_NAME','?.')
4577	> IF ( ok_var(var_name) ) THEN
4578	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
4579	> IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
4580	> zdz1(:,:)=temp1(:,:)
4581	> ELSE
4582	> i_notfound = i_notfound + 1
4583	> ENDIF
4584	> ENDIF
4585	>
4586	> var_name= 'zdz2'
4587	> CALL ioconf_setatt('UNITS', '-')
4588	> CALL ioconf_setatt('LONG_NAME','?.')
4589	> IF ( ok_var(var_name) ) THEN
4590	> CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
4591	> IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
4592	> zdz2(:,:)=temp(:,:)
4593	> ELSE
4594	> i_notfound = i_notfound + 1
4595	> ENDIF
4596	> ENDIF
4597	>
4598	> var_name='temp_sol_beg'
4599	> CALL ioconf_setatt('UNITS', 'K')
4600	> CALL ioconf_setatt('LONG_NAME','Old Surface temperature')
4601	> IF ( ok_var(var_name) ) THEN
4602	> CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
4603	> IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
4604	> temp_sol_beg(:) = temp2(:)
4605	> ELSE
4606	> i_notfound = i_notfound + 1
4607	> ENDIF
4608	> ENDIF
4609	>
4610	> !
4611	477a500,508
4612	> !
4613	> ! Should the restart be incomplete for some reason we reset ptn as well
4614	> !
4615	> IF ( i_notfound > 0 ) THEN
4616	> ptn_default = 280._r_std
4617	> CALL getin_p('THERMOSOIL_TPRO',ptn_default)
4618	> ptn(:,:) = ptn_default
4619	> ENDIF
4620	> !!$ ptn(:,:) = 234.900240607854926
4621	Seulement dans /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba: thermosoil.f90.old
4622	_______________________________________________________________________________________________________________________
4623	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE_OL/dim2_driver.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE_OL/dim2_driver.f90
4624	369c369
4625	< & 'Restart file initialized : itau_dep, itau_fin, date0_rest, dt_rest', &
4626	---
4627	> & 'Restart file initialized : itau_dep, itau_fin, date0_rest, dt_rest : ', &
4628	542a543,546
4629	> !
4630	> ! Jan : Quick fix until the steping is solved correctly in this version :-(
4631	> !
4632	> for_offset = -itau_dep
4633	586c590
4634	< STOP 'dim2_driver'
4635	---
4636	> !! STOP 'dim2_driver'
4637	886a891
4638	> !
4639	1090a1096,1100
4640	> ! petAcoef(:,:) = zero
4641	> ! petBcoef(:,:) = zero
4642	> ! peqAcoef(:,:) = zero
4643	> ! peqBcoef(:,:) = zero
4644	>
4645	1197a1208
4646	> !
4647	1657c1668,1669
4648	< CALL Write_Load_Balance(Get_cpu_time(timer_mpi))
4649	---
4650	> !! CALL Write_Load_Balance(Get_cpu_time(timer_mpi))
4651	> CALL Write_Load_Balance(REAL(Get_cpu_time(timer_mpi), r_std))
4652	_______________________________________________________________________________________________________________________
4653	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE_OL/readdim2.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE_OL/readdim2.f90
4654	68c68
4655	< LOGICAL :: debug = .FALSE.
4656	---
4657	> LOGICAL :: debug = .TRUE.
4658	110c110,111
4659	< ! WRITE(numout,*) "forcing_info : Forcing time step out of flinopen : ",dt_force
4660	---
4661	> IF ( debug ) WRITE(numout,*) "forcing_info : first and last itau : ", itau(1), itau(ttm_full)
4662	> IF ( debug ) WRITE(numout,*) "forcing_info : Forcing time step out of flinopen : ",dt_force
4663	485c486
4664	< & fcontfrac
4665	---
4666	> & fcontfrac, tmp_cont
4667	493a495
4668	> INTEGER,DIMENSION(2) :: kloc
4669	557a560,561
4670	> tmp_cont(:,:) = zero
4671	> !-
4672	564a569,600
4673	> DO ik=1,nbindex
4674	> j = ((kindex(ik)-1)/iim) + 1
4675	> i = (kindex(ik) - (j-1)*iim)
4676	> tmp_cont(i,j) = un
4677	> ENDDO
4678	> !-
4679	> !- Intercept some stupid forcing variables
4680	> !-
4681	> IF ( MINVAL(precip) < zero ) THEN
4682	> kloc=MINLOC(precip)
4683	> WRITE(,) "Negative precipitation in forcing_read :", kloc(1), kloc(2), precip(kloc(1), kloc(2))
4684	> STOP
4685	> ENDIF
4686	> !-
4687	> IF ( MINVAL(snowf) < zero ) THEN
4688	> kloc=MINLOC(snowf)
4689	> WRITE(,) "Negative precipitation in forcing_read :", kloc(1), kloc(2), snowf(kloc(1), kloc(2))
4690	> STOP
4691	> ENDIF
4692	> !-
4693	> IF ( MINVAL(pb) < 10000 ) THEN
4694	> kloc=MINLOC(pb)
4695	> WRITE(,) "Weak surface pressure (< 100hPa) in forcing_read :", kloc(1), kloc(2), pb(kloc(1), kloc(2))
4696	> STOP
4697	> ENDIF
4698	> !-
4699	> IF ( ABS(MAXVAL(utmp_cont)) > 100 .OR. ABS(MAXVAL(vtmp_cont)) > 100 ) THEN
4700	> WRITE(numout,*) "Too strong wind at the surface (u or v > 100 m/s) :", &
4701	> MAXVAL(utmp_cont), MAXVAL(vtmp_cont)
4702	> STOP
4703	> ENDIF
4704	> !
4705	678,679c714
4706	< ! to bypass FPE problem on integer convertion with missing_value heigher than precision
4707	< INTEGER, PARAMETER :: undef_int = 999999999
4708	---
4709	> !
4710	720c755
4711	< IF (is_watchout) THEN
4712	---
4713	> !! IF (is_watchout) THEN
4714	727c762
4715	< ENDIF
4716	---
4717	> !! ENDIF
4718	950c985
4719	< & 'The dates : ',itau_read,itau_split,itau_read_nm1,itau_read_n
4720	---
4721	> & 'The dates : ',itau_read,last_read,itau_split,itau_read_nm1,itau_read_n
4722	963a999
4723	> IF (check) WRITE(numout,*) "Date nm1 : ", itau_read_nm1, itau_read
4724	1356c1392
4725	< IF (MINVAL(tair(:,:)) < 100.) THEN
4726	---
4727	> IF (MINVAL(tair(:,:)) < 100. .AND. MINVAL(tair(:,:)) < 0.) THEN
4728	1381c1417
4729	< IF (tair(i,j) < 500.) THEN
4730	---
4731	> IF (tair(i,j) < 500. .AND. tair(i,j) > 5.) THEN
4732	_______________________________________________________________________________________________________________________
4733	diff -w --ignore-all-space --ignore-case --recursive --exclude='cvs_diff' --exclude='.flc' --exclude='.bak' --exclude='.svn' --exclude='.lst' --exclude='i..L' --exclude='~' --exclude=Entries --exclude=Tag --exclude=CVS --exclude Makefile /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE_OL/weather.f90 /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE_OL/weather.f90
4734	1490,1491c1490,1491
4735	< INTEGER,DIMENSION(iim),INTENT(OUT) :: iind
4736	< INTEGER,DIMENSION(jjm),INTENT(OUT) :: jind
4737	---
4738	> INTEGER,DIMENSION(iim),INTENT(INOUT) :: iind
4739	> INTEGER,DIMENSION(jjm),INTENT(INOUT) :: jind
4740	1671c1671
4741	< WRITE(numout,'(a1,$)') maskchar(NINT(lsm_file(i,j)))
4742	---
4743	> WRITE(numout,*) maskchar(NINT(lsm_file(i,j)))
4744	1730c1730
4745	< WRITE(numout,'(a1,$)') maskchar(NINT(mask(i,j)))
4746	---
4747	> WRITE(numout,*) maskchar(NINT(mask(i,j)))
4748	3058a3059
4749	> REAL,DIMENSION(1) :: tau_tmp
4750	3063c3064,3065
4751	< iret = NF90_PUT_VAR (dump_id,timestp_id,(/ REAL(itau) /), &
4752	---
4753	> tau_tmp = REAL(itau)
4754	> iret = NF90_PUT_VAR (dump_id,timestp_id, tau_tmp, &
4755	3065c3067,3068
4756	< iret = NF90_PUT_VAR (dump_id,time_id,(/ REAL(itau)*dt_force /), &
4757	---
4758	> tau_tmp = REAL(itau)*dt_force
4759	> iret = NF90_PUT_VAR (dump_id,time_id, tau_tmp, &

Download in other formats:

Original Format