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, 10 years ago)

Fichier des diffs entre la version LMD et la version 1.9.4.1

Line 
1Seulement dans /login/IPSL_CODE/ORCHIDEE_CVS_1_9_4_1/modeles/ORCHIDEE: Doxyfile_ORCHIDEE
2_______________________________________________________________________________________________________________________
3diff -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
42c2
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 $
840a41
9> #-Q- sx9mercure       mv $*.mod $(MODDIR)
10_______________________________________________________________________________________________________________________
11diff -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
1226,27d25
13<   ! earth radius
14<   REAL(r_std), PARAMETER :: R_Earth = 6378000.
15_______________________________________________________________________________________________________________________
16diff -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
1734,35d33
18<   REAL(r_std), PARAMETER                          :: R_Earth = 6378000.
19<   REAL(r_std), PARAMETER                          :: mincos  = 0.0001
20339a338,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
26357c360
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)
30383a387,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>
43440c455
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
47442c457,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)
54444c462,463
55<
56---
57>     INTEGER        :: ALLOC_ERR
58>     !
59472a492,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_minlon*incp)-2:NINT(domain_maxlon*incp)+2, &
92>          & NINT(domain_minlat*incp)-2:NINT(domain_maxlat*incp)+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_low*incp):NINT(lon_up*incp),NINT(lat_low*incp):NINT(lat_up*incp))=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>     !
120491,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.)
126494,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.)
132505c584
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'
136508,511d586
137< #ifdef INTERPOL_ADVANCE
138<     WRITE(numout,'(2a40)')'0%--------------------------------------', &
139<          & '------------------------------------100%'
140< #endif
141518,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)
170547,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.
175553a597
176>        !
177562,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)
183567,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)
189571,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(lonrel*incp),NINT(latrel*incp))
199577,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
230603a623
231>        pp = fine_ind(NINT(louprel*incp),NINT(lauprel*incp))
232604a625,627
233>        IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
234>           pos_pnt(ipp+1) = pp
235>           ipp = ipp + 1
236606,610c629
237<                 ENDIF       ! IF lat
238<              ENDIF          ! IF lon
239<           ENDDO
240<           !
241<           ! De-select the marked points
242---
243>        pp = fine_ind(NINT(louprel*incp),NINT(lalowrel*incp))
244612,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
255619c635
256<           ENDDO
257---
258>        pp = fine_ind(NINT(lolowrel*incp),NINT(lauprel*incp))
259621,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
277634a641
278>        pp = fine_ind(NINT(lolowrel*incp),NINT(lalowrel*incp))
279636,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
297650,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)
310663,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)
315667a649
316>        IF ( ipp > zero ) THEN
317669,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)
327689a668,669
328>              !
329>              !
330701a682
331>                    !
332704d684
333<                       EXIT loopnbind
334715d694
335<                       found = .TRUE.
336718,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
347737,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
352742,743d717
353<        ENDIF
354<     ENDIF
355745c719,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
363749a728
364>     DEALLOCATE (fine_ind)
365817a797,798
366>     INTEGER(i_std) :: ib, iv
367>
368820a802,805
369>
370>     sub_index(:,:,:) = 0
371>     sub_area(:,:) = zero
372>
373825d809
374<   
375_______________________________________________________________________________________________________________________
376diff -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
3772c2
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 $
38140a41
382> #-Q- sx9mercure       mv $*.mod $(MODDIR)
383_______________________________________________________________________________________________________________________
384diff -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
385128a129,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_______________________________________________________________________________________________________________________
404diff -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
405569c569
406<     LOGICAL, PARAMETER :: check=.FALSE.
407---
408>     LOGICAL, PARAMETER :: check=.TRUE.
409_______________________________________________________________________________________________________________________
410diff -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
4112c2
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 $
41542a43
416> #-Q- sx9mercure       mv $*.mod $(MODDIR)
417_______________________________________________________________________________________________________________________
418diff -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
41937a38,40
420>   REAL(r_std), PARAMETER                          :: R_Earth = 6378000.
421>   REAL(r_std), PARAMETER                          :: mincos  = 0.0001
422> !
42344a48
424>   INTEGER(i_std),PARAMETER :: undef_int = 999999999
42571a76,77
426>     LOGICAL :: do_floodplains
427>     LOGICAL :: do_irrigation
428_______________________________________________________________________________________________________________________
429diff -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
43011a12,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> !
43523a28,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
44825a42,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
45590a113,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
47092c128
471<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: nvan = &
472---
473>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: nvan_fao = &
47494,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 = &
484101,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 /)
494108c136
495<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcr = &
496---
497>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcr_fao = &
498111c139
499<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcs = &
500---
501>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcs_fao = &
502113,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
507118c142
508<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: ks = &
509---
510>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: ks_fao = &
511120,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 = &
517124c148
518<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: free_drain_max = &
519---
520>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: free_drain_max_fao = &
521127c151
522<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcf = &
523---
524>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcf_fao = &
525130c154
526<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mcw = &
527---
528>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mcw_fao = &
529131a156
530> ! & (/ 0.07_r_std, 0.085_r_std, 0.10_r_std /)
531133c158
532<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mc_awet = &
533---
534>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mc_awet_fao = &
535136c161
536<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: mc_adry = &
537---
538>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: mc_adry_fao = &
539139c164
540<   REAL(r_std),PARAMETER,DIMENSION(nstm) :: psis = &
541---
542>   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: psis_fao = &
543140a166,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_______________________________________________________________________________________________________________________
628diff -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
629173,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.
635179,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', &
641233c228
642<   INTEGER(i_std),DIMENSION(nvm,nstm) ::  pref_soil_veg
643---
644>   INTEGER(i_std),DIMENSION(nvm) ::  pref_soil_veg
645_______________________________________________________________________________________________________________________
646diff -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
6472c2
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 $
65184a85
652> #-Q- sx9mercure       mv $*.mod $(MODDIR)
653_______________________________________________________________________________________________________________________
654diff -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
6551d0
656<
65738c37
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
661185a185,187
662>     LOGICAL, SAVE                                         :: fatmco2       !! Flag to force the value of atmospheric CO2 for vegetation.
663>     REAL(r_std), SAVE                                     :: atmco2        !! atmospheric CO2
664>     !
665225a228,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>        !
693287d315
694<        zccanopy(ik)     = ccanopy(i,j)
695296a325,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>     !
707362c401
708<     IF ( check ) WRITE(numout,*) 'Calling sechiba'
709---
710>     IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_main_2d'
711364a404
712>        !
713370c410
714<        CALL histsync(hist_id)
715---
716> !!$       CALL histsync(hist_id)
717375c415
718<        CALL histsync(hist2_id)
719---
720> !!$       CALL histsync(hist2_id)
721378a419
722>     !
723383,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, &
729467,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>           !
735486,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)
738495a531,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)
741506,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
749640a678,680
750>     LOGICAL, SAVE                                         :: fatmco2       !! Flag to force the value of atmospheric CO2 for vegetation.
751>     REAL(r_std), SAVE                                     :: atmco2        !! atmospheric CO2
752>     !
753695a736,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>        !
780733d798
781<        zccanopy(ik)     = ccanopy(kindex(ik))
782742a808,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
791747c820
792<     IF ( check ) WRITE(numout,*) 'Calling sechiba'
793---
794>     IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_main_1d'
795755c828
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)
799760c833
800<        CALL histsync(hist_id)
801---
802> !!$       CALL histsync(hist_id)
803763c836
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)
807765c838
808<        CALL histsync(hist2_id)
809---
810> !!$       CALL histsync(hist2_id)
811773,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, &
817856,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>           !
823875,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)
826884a952,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)
8291233d1301
830<       
8311323c1391
832<     IF ( check ) WRITE(numout,*) 'Calling sechiba'
833---
834>     IF ( check ) WRITE(numout,*) 'Calling sechiba intersurf_gathered'
8351370c1438
836<        CALL histsync(hist_id)
837---
838> !!$       CALL histsync(hist_id)
8391375c1443
840<        CALL histsync(hist2_id)
841---
842> !!$       CALL histsync(hist2_id)
8431383,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, &
8491486,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)
8521505,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)
8551515a1578,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)
8581675c1739,1741
859<
860---
861>     ! LOCAL
862>     CHARACTER(LEN=30) :: classif
863>     !
8641685a1752,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>     !
8891743,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>     !
9331890a2011
934>     !
9352031c2152
936<     INTEGER(i_std)                         :: ier
937---
938>     INTEGER(i_std)                         :: ier, jv
9392104a2226,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>     !
9452241c2367
946<                &    nslm, solay, solayax_id)
947---
948>                &    nslm, diaglev(1:nslm), solayax_id)
9492261a2388,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>        !
9622289,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)
9782356a2491,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
9832369a2508,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>             
9892371a2515,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)
9932386a2533,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)
9982443a2594,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
10072445c2603
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)
10112447,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>           !
10292464c2632
1030< !MM
1031---
1032>        !
10332470a2639,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)
10362473,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)
10392476a2645,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
10512528c2707
1052<                &    nslm, solay, solayax_id)
1053---
1054>                &    nslm, diaglev(1:nslm), solayax_id)
10552582c2761
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',&
10592585a2765,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
10682590a2778,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)
10712595c2784,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', &
10792615a2809,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)
10822621a2817,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)
10852627a2825,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)
10972641,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)
11032645,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     
11172663a2880,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)
11202667a2886,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)
11272971a3196,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>           !
11322986a3215,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)
11373004a3237,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)
11403041a3276,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)
11473045a3286,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)
11503049a3292,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)
11553070c3316,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', '%', &
11613071a3320,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)
11683124c3378
1169<                   &    nslm, solay, solayax_id2)
1170---
1171>                   &    nslm, diaglev(1:nslm), solayax_id2)
11723130c3384
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)
11763132c3386
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)
11803134c3388
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)
11843153c3407
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)
11883155c3409
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)
11923157c3411
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)
11963159c3413
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)
12003176,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)
12083180c3434,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)
12143187c3443
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)
12183188a3445,3446
1219>                & iim,jjm, hori_id2, 1,1,1, -99, 32, avescatter2(7), dt, dw2)
1220>           CALL histdef(hist2_id, 'SWI', 'Soil wetness index','-',  &
12213190c3448
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',  &
12253191a3450,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)
12283197c3457
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)
12323200c3460
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)
12363205c3465
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)
12403212c3472
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)
12443214c3474
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)
12483216c3476,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)
12543222c3484
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)
12583228c3490
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)
12623230c3492
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)
12663232c3494
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)
12703237c3499,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)
12783239,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)
12843247c3511
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)
12883251c3515
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)
12923258a3523,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)
12963368a3636
1297>
1298_______________________________________________________________________________________________________________________
1299diff -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
13004c4
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 $
130431,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<
131650c39,40
1317<
1318---
1319>   REAL(r_std), SAVE                                :: slope_default = 0.1
1320>   !
132156d45
1322<   LOGICAL, SAVE                                   :: old_lai = .FALSE.         ! Old Lai Map interpolation
132358d46
1324<   LOGICAL, SAVE                                   :: old_veget = .FALSE.         ! Old veget Map interpolation
132561a50
1326>   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilclass_default
132771c60
1328<        IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
1329---
1330>        IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
133176c65
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, &
1335117a107
1336>     INTEGER(i_std), DIMENSION(kjpindex), INTENT(inout)  :: njsc             !! Soil type map to be created from the input map
1337118a109
1338>     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: frac_bare        !! Bare soil fraction for each tile
1339124c115,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
1344157,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>
1351244a238,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>        !
1355248,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)
1367294c296
1368<     IF (check) WRITE(*,*) "slowproc: day_counter 3",day_counter
1369---
1370>     IF (check) WRITE(numout,*) "slowproc: day_counter 3",day_counter
1371308c310
1372<        IF (check) WRITE(*,*) "slowproc: day_counter 2",day_counter
1373---
1374>        IF (check) WRITE(numout,*) "slowproc: day_counter 2",day_counter
1375360a363,365
1376>     IF (EndOfYear .AND. .NOT. land_use_updated) THEN
1377>        lcchange=.FALSE.
1378>     ENDIF
1379410c415
1380<                lalo,resolution,lai,month,day,read_lai,laimap)
1381---
1382>                lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
1383450,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)
1389472a478
1390>     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)    :: frac_bare      !! Bare soil fraction for each tile
1391478c484,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
1397486c494,495
1398<     INTEGER(i_std)                                       :: ji, jv, ier
1399---
1400>     INTEGER(i_std)                                       :: ji, jv, ier, jst
1401> !
1402487a497,498
1403>     LOGICAL                                              :: get_slope
1404> !
1405490a502,503
1406>     REAL(r_std), DIMENSION(kjpindex)                     :: tmp_real
1407>
1408491a505
1409>     REAL(r_std), DIMENSION (kjpindex,nscm)               :: soilclass   !! fraction of soil type
1410502a517,518
1411>     IF (long_print) WRITE (numout,*) "In slowproc_init"
1412>
1413509a526,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>     !
1421537a561,563
1422>     !
1423>     IF (long_print) WRITE (numout,*) 'slowproc_init : End of allocations'
1424>
1425556a583,584
1426>     IF (long_print) WRITE (numout,*) 'slowproc_init : End of allocationsGot restart file set'
1427>
1428612c640
1429<        CALL getin('VEGET_YEAR', veget_year_orig)
1430---
1431>        CALL getin_p('VEGET_YEAR', veget_year_orig)
1432655c683,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'
1438658c688,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)
1447664a700,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>     !
1458669a715,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>     !
1464703,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>     !
1484782c818,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
1496853c897
1497<        CALL setvar_p (soiltype, val_exp, 'SOIL_FRACTIONS', soiltype_default)
1498---
1499>        CALL setvar_p (soilclass, val_exp, 'SOIL_FRACTIONS', soilclass_default)
1500862a907,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>
1510874a928,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
1526884,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<
1538897,898c954
1539<                 IF ( .NOT. old_veget ) THEN
1540<                    ! NEW slowproc interpol :
1541---
1542>                 ! slowproc interpol :
1543900,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
1548986,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<
15601002,1003d1042
1561<              IF ( .NOT. old_lai ) THEN
1562<                 ! NEW slowproc interlai :
15631005,1008d1043
1564<              ELSE
1565<                 ! OLD slowproc interlai :
1566<                 CALL slowproc_interlai(kjpindex, lalo,  resolution, laimap)
1567<              ENDIF
15681014c1049
1569<                   lalo,resolution,lai,month,day,read_lai,laimap)
1570---
1571>                   lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
15721030,1031d1064
1573<                 IF ( .NOT. old_lai ) THEN
1574<                    ! NEW slowproc interlai :
15751033,1036c1066
1576<                 ELSE
1577<                    ! OLD slowproc interlai :
1578<                    CALL slowproc_interlai(kjpindex, lalo,  resolution, laimap)
1579<                 ENDIF
1580---
1581>                 !
15821042c1072
1583<                lalo,resolution,lai,month,day,read_lai,laimap)
1584---
1585>                lalo, resolution, lai, frac_bare, month, day, read_lai, laimap)
15861061c1091,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.&
15981062a1101
1599>           CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
16001064c1103,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)
16421069,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<     !
16711103a1153,1155
1672>     !
1673>     ! Clean-up
1674>     !
16751122a1175,1178
1676>     IF ( ALLOCATED (soilclass_default)) DEALLOCATE (soilclass_default)
1677>     !
1678>
1679>     !
16801489c1545
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)
16841505a1562
1685>     REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)   :: frac_bare  !! Bare fraction for each tile
16861508a1566
1687>     REAL(r_std), SAVE                                  :: extcoef     !! Extinction coefficient for bare soil with LAI
16881509a1568,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
16991511,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
17071523d1586
1708<           IF ( .NOT. read_lai ) THEN
17091525,1529d1587
1710<           ELSE
1711<              DO ji = 1,kjpindex
1712<                 lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
1713<              ENDDO
1714<           ENDIF
17151535d1592
1716<           IF ( .NOT. read_lai ) THEN
17171538a1596,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>        !
17311539a1610,1614
1732>        !
1733>        DO jv = 1,nvm
1734>           !
1735>           ! If January
1736>           !
17371545a1621,1623
1738>           !
1739>           ! If December
1740>           !
17411551a1630,1632
1742>           !
1743>           ! ELSE
1744>           !
17451559,1563c1640
1746<           ENDIF
1747<           !
1748<        CASE default
1749<           !
1750<           ! 3. Problem
1751---
1752>        ENDDO
17531565,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
17611570a1644,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))
17671571a1650
1768>     ENDIF
17691574a1654
1770>   !!
17711576d1655
1772< !MM TAG 1.6 model !
17731578c1657
1774<   SUBROUTINE slowproc_interlai_OLD(nbpt, lalo,  resolution, laimap)
1775---
1776>   SUBROUTINE slowproc_interlai(nbpt, lalo,  resolution, neighbours, contfrac, laimap)
17771587a1667,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>     !
17811594,1595d1675
1782<     REAL(r_std), PARAMETER                          :: R_Earth = 6378000., min_sechiba=1.E-8
1783<     !
17841598,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
17981606,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
18181617c1693
1819<     pi = 4. * ATAN(1.)
1820---
1821>     INTEGER                  :: ALLOC_ERR
18221621c1697
1823<     !Config If   = !LAI_MAP
1824---
1825>     !Config If   = LAI_MAP
18261637,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
18561642d1737
1857<     WRITE(numout,*) 'slowproc_interlai : Reading the LAI file'
18581648c1743,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
18651655,1657d1752
1866<     CALL bcast(mask_lu)
1867<     !
1868<     WRITE(numout,*) 'slowproc_interlai : ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml)
18691658a1754,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
18821660,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(:)
18991673a1770,1771
1900>     DO jp=1,jml
1901>        lon(:,jp) = lon_lu(:)
19021676,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
19141681,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
19211686c1783
1922<     DO ib =1, nbpt
1923---
1924>     mask(:,:) = 0
19251688,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
19321692c1788,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.
19371694,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
19441697c1795
1945<        coslat = pi/180. * R_Earth
1946---
1947>           ! Resolution in longitude
19481699,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
19601701a1806
1961>           ! Resolution in latitude
19621702a1808,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
19701705,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
20081741c1819,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.
20131743,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
20281752a1827
2029>     CALL bcast(nbvmax)
20301754,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'
20501772c1831
2051<              ! We were lucky
2052---
2053>     ok_interpol = .FALSE.
20541774,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
23372072d1853
2338<       
23392074,2080d1854
2340<        IF ( .NOT. ok_interpol ) THEN
2341<           DEALLOCATE(sub_area)
2342<           DEALLOCATE(sub_index)
2343<        ENDIF
2344<        !
2345<        nbvmax = nbvmax * 2
2346<     ENDDO
23472102,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
23572121a1887,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>     !
23742133a1915
2375>     DEALLOCATE (resol_lu)
23762137c1919
2377<   END SUBROUTINE slowproc_interlai_NEW
2378---
2379>   END SUBROUTINE slowproc_interlai
23802546c2328
2381<              WRITE(*,*) "SUM(veget_last(",ib,")) = ",sum_veg
2382---
2383>              WRITE(numout,*) "SUM(veget_last(",ib,")) = ",sum_veg
23842718d2499
2385< !MM TAG 1.6 model !
23862720c2501
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 )
23902730a2512
2391>     REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac         !! Fraction of continent in the grid
23922736a2519,2520
2393>     LOGICAL ::           ok_interpol  ! optionnal return of aggregate_vec
2394>     !
23952739d2522
2396<     REAL(r_std), PARAMETER                          :: R_Earth = 6378000., min_sechiba=1.E-8
23972744,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
24032748,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
24122754,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
24242760c2542,2543
2425<     pi = 4. * ATAN(1.)
2426---
2427>     n_origveg(:,:) = zero
2428>     n_found(:) = zero
24292766a2550
2430>     !Config If   = !LAND_USE
24312783,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)
24682803,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'
24762807,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
24822812,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.
25392865c2615,2616
2540<     DO ip=1,iml
2541---
2542>     resol_lon = 5000.
2543>     resol_lat = 5000.
25442867c2618,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.
25492869,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
25642878a2626
2565>     CALL bcast(nbvmax)
25662880,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"
26062918c2630,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
26132920,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
26242924a2641,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
26332926,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)
26462936,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)
26532939d2656
2654<             ENDDO
26552940a2658,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)
26622942d2664
2663<             !
26642944,2947d2665
2665<           !
2666<         ENDIF ! lucky/not lucky
2667<         !
2668<       ENDIF     ! in the domain
26692949c2667
2670<
2671---
2672>     !
26732954,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<
26832968c2677
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(:)
26872973d2681
2688<
26892977a2686,2688
2690>     veget(:,:) = zero
2691>     frac_nobio(:,:) = zero
2692>     !
26932990,2992c2701
2694<     !
2695<     WRITE(numout,*)
2696<     WRITE(numout,*) 'Interpolation Done'
2697---
2698>     WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
26993024,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)
27073062,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)
27143068a2772,2775
2715>     DEALLOCATE(lat_ful, lon_ful)
2716>     DEALLOCATE(sub_index)
2717>     DEALLOCATE(sub_area)
2718>
27193072c2779
2720<  END SUBROUTINE slowproc_interpol_OLD
2721---
2722>   END SUBROUTINE slowproc_interpol
27233076c2783
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 )
27273113c2820
2728<     INTEGER(i_std) :: idi, jv, inear, nbvmax
2729---
2730>     INTEGER(i_std) :: idi, jv, inear, nbvmax, nix, njx
27313133c2840
2732<     CALL getin_p('VEGETATION_FILE',filename)
2733---
2734>     CALL getin('VEGETATION_FILE',filename)
27353135,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)
27433163d2865
2744<     IF (is_root_prc) THEN
27453167a2870,2872
2746>     WRITE(numout,*) 'File name : ', filename
2747>     WRITE(numout,*) 'Min and max vegetation numbers : ', MINVAL(vegmap), MAXVAL(vegmap)
2748>     !
27493169,3174d2873
2750<     ENDIF
2751<     
2752<     CALL bcast(lon_ful)
2753<     CALL bcast(lat_ful)
2754<     CALL bcast(vegmap)
2755<     
27563193,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
27703198a2904,2905
2771>     ! We are on a mono proc routine and thus nbvmax does not need to be broadcasted
2772>     !
27733200d2906
2774<     DO WHILE ( .NOT. ok_interpol )
27753202c2908
2776<        WRITE(numout,*) "nbvmax = ",nbvmax
2777---
2778>     WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
27793219c2925
2780<        CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, &
2781---
2782>     CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, &
27833223,3228d2928
2784<        IF ( .NOT. ok_interpol ) THEN
2785<           DEALLOCATE(sub_area)
2786<           DEALLOCATE(sub_index)
2787<           !
2788<           nbvmax = nbvmax * 2
2789<        ELSE
27903240,3241d2939
2791<        ENDIF
2792<     ENDDO
27933353c3051,3052
2794<   END SUBROUTINE slowproc_interpol_NEW
2795---
2796>   END SUBROUTINE slowproc_interpol_g
2797>
27983356,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( cospa*coscolat + sinpa*sincolat*COS(p))
2842>        cosang(i) = cospa*coscolat + sinpa*sincolat*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
28543359c3108
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)
28583361a3111,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>     !
28643369a3124
2865>     REAL(r_std), INTENT(in)       :: contfrac(nbpt)     ! Fraction of land in each grid box.
28663373,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
28723376d3130
2873<     !  0.3 LOCAL
28743378,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
28793380a3134
2880>     INTEGER(i_std)               :: nbvmax
28813383,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
28873386,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
29123399c3152
2913<     pi = 4. * ATAN(1.)
2914---
2915>     ! Number of texture classes in Zobler
29163401c3154,3155
2917<     CALL get_vegcorr (nolson,vegcorr,nobiocorr)
2918---
2919>     INTEGER(i_std), PARAMETER :: nzobler = 7
2920>     REAL(r_std),ALLOCATABLE   :: textfrac_table(:,:)
29213403,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
29323412,3413c3159
2933<     filename = '../surfmap/carteveg5km.nc'
2934<     CALL getin('VEGETATION_FILE',filename)
2935---
2936>     INTEGER                  :: ALLOC_ERR
29373415d3160
2938<     CALL flininfo(filename, iml, jml, lml, tml, fid)
29393416a3162
2940>     !  Needs to be a configurable variable
29413418,3420d3163
2942<     ALLOCATE(lat_ful(iml))
2943<     ALLOCATE(lon_ful(iml))
2944<     ALLOCATE(vegmap(iml))
29453422c3165,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.
29573424,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)
29643428,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)
29753431,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)
29923439a3194,3197
2993>     CALL bcast(iml)
2994>     CALL bcast(jml)
2995>     CALL bcast(lml)
2996>     CALL bcast(tml)
29973441,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.
30043446c3201,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
30433448,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>     
30553452c3246,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)
30603454,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
30683457c3254
3069<        coslat = pi/180. * R_Earth
3070---
3071>     IF (is_root_prc) CALL flinclo(fid)
30723459,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
30773463,3464c3259
3078<        veget(ib,:) = 0.0
3079<        frac_nobio (ib,:) = 0.0
3080---
3081>     ! Mask of permitted variables.
30823466c3261,3263
3083<     ENDDO
3084---
3085>     mask(:,:) = zero
3086>     DO ip=1,iml
3087>        DO jp=1,jml
30883468c3265,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
30943470,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
31093483,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
31223487,3489c3280
3123<     WRITE(numout,*) 'Interpolating the vegetation map :'
3124<     WRITE(numout,'(2a40)')'0%--------------------------------------', &
3125<                    & '------------------------------------100%'
3126---
3127>           ! Resolution in latitude
31283491,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
31393494c3290,3291
3140<     DO ip=1,iml
3141---
3142>        ENDDO
3143>     ENDDO
31443496c3293,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.
31493498,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
31643507a3301
3165>     CALL bcast(nbvmax)
31663509,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"
31713512,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.
31783517,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
31863521a3310,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
32133523,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)
32203526,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
32253531c3342
3226<           ! We were lucky
3227---
3228>     ! Depending on the selected classification, the interpolation will be different
32293533c3344,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')
32403535c3353
3241<         ELSE
3242---
3243>        soilclass_default=soilclass_default_fao
32443537,3538c3355
3245<           ! Otherwise, look everywhere.
3246<           ! Begin close to last grid point.
3247---
3248>        PRINT *, "Using a soilclass map with Zobler classification"
32493540,3541c3357
3250<           found = .FALSE.
3251<           idi = 1
3252---
3253>        ALLOCATE(textfrac_table(nzobler,ntext))
32543543c3359
3255<           DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
3256---
3257>        CALL get_soilcorr (nzobler, textfrac_table)
32583545d3360
3259<             ! forward and backward
32603547c3362
3261<             DO ii = 1,2
3262---
3263>        DO ib =1, nbpt
32643549,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
32723555,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
32783561,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)
32843565,3566c3369
3285<                 ENDIF
3286<               ENDIF
3287---
3288>           !    Check that we found some points
32893568c3371,3372
3290<             ENDDO
3291---
3292>           soilclass(ib,:) = 0.0
3293>           clayfraction(ib) = 0.0
32943570c3374,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
33023571a3380,3381
3303>              DO ilf = 1,fopt
3304>                 solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
33053574c3384
3306<         ENDIF ! lucky/not lucky
3307---
3308>              sgn = zero
33093576,3578c3386
3310<       ENDIF     ! in the domain
3311<     ENDDO
3312<
3313---
3314>              !   Compute the average bare soil albedo parameters
33153580,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
33203583c3390
3321<
3322---
3323>                 ! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean
33243585,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
33563588,3589d3420
3357<     DO ib = 1, nbpt
3358<       n_found(ib) = SUM( n_origveg(ib,:) )
33593591,3593d3421
3360<
3361<     !
3362<     ! determine fraction of Olson vegetation type in each box of the coarse grid
33633595,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
33743604,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
33863607c3434
3387<     DO vid = 1, nolson
3388---
3389>           ENDIF
33903609,3610d3435
3391<       DO jv = 1, nvm
3392<         veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
33933611a3437
3394>   
33953613,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")
34013617c3441
3402<     ENDDO
3403---
3404>        soilclass_default=soilclass_default_fao
34053618a3443
3406>        PRINT *, "Using a soilclass map with fao classification"
34073620,3621c3445
3408<     WRITE(numout,*)
3409<     WRITE(numout,*) 'Interpolation Done'
3410---
3411>        ALLOCATE(textfrac_table(nscm,ntext))
34123623c3447
3413<     !   Clean up the point of the map
3414---
3415>        CALL get_soilcorr (nscm, textfrac_table)
34163627,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
34223631d3452
3423<           ! Now we need to handle some exceptions
34243633,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)
34323639,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
34403645,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
34493650a3461,3464
3450>           IF ( fopt .EQ. 0) THEN
3451>              nbexp = nbexp + 1
3452>              soilclass(ib,:) = soilclass_default(:)
3453>              clayfraction(ib) = clayfraction_default
34543653,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
34623657,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)
34683663,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
34743667,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
34823671d3477
3483<           ENDIF
34843672a3479,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
34953674a3491
3496>              ENDDO
34973676c3493
3498<        !  Limit the smalest vegetation fraction to 0.5%
3499---
3500>              ! Normalize the surface
35013678,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
35133682,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
35193687a3504
3520>           ENDIF
35213691,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)
35293699c3509
3530<     RETURN
3531---
3532>     CASE("fao2")
35333701,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
35413706a3513
3542>        PRINT *, "Using a soilclass map with fao2 classification"
35433707a3515
3544>        ALLOCATE(textfrac_table(nscm,ntext))
35453709c3517
3546<     !  0.1 INPUT
3547---
3548>        CALL get_soilcorr (nscm, textfrac_table)
35493711,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
35693728a3521
3570>           ! GO through the point we have found
35713730,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
35853744c3524
3586<     INTEGER                  :: ALLOC_ERR
3587---
3588>           fopt = COUNT(sub_area(ib,:) > zero)
35893746,3747c3526
3590<     n_origveg(:,:) = zero
3591<     n_found(:) = zero
3592---
3593>           !    Check that we found some points
35943749c3528,3529
3595<     CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3596---
3597>           soilclass(ib,:) = 0.0
3598>           clayfraction(ib) = 0.0
35993751,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
36153761,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
36233764d3541
3624<     CALL flininfo(filename, iml, jml, lml, tml, fid)
36253765a3543
3626>              !   Compute the average bare soil albedo parameters
36273767,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
36483786c3547
3649<     WRITE(numout,*) 'Reading the OLSON type vegetation file'
3650---
3651>              DO ilf = 1,fopt
36523788,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)
36563792d3549
3657<     CALL flinclo(fid)
36583794,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,2*solt(ilf)) = soilclass(ib,2*solt(ilf)) + sub_area(ib,ilf)
3671>                    ELSE
3672>                       soilclass(ib,2*solt(ilf)-1) = soilclass(ib,2*solt(ilf)-1) + sub_area(ib,ilf)
36733803,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'
37053828,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
37123835d3564
3713<        sub_area(:,:)=zero
37143837,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
37203841,3843c3568
3721<        IF ( .NOT. ok_interpol ) THEN
3722<           DEALLOCATE(sub_area)
3723<           DEALLOCATE(sub_index)
3724---
3725>              ! Normalize the surface
37263845c3570,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
37333847,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
37483859,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<     !
37543865c3579
3755<     ! determine fraction of Olson vegetation type in each box of the coarse grid
3756---
3757>           ENDIF
37583867,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
37653875,3876d3582
3766<     ! now finally calculate coarse vegetation map
3767<     ! Find which model vegetation corresponds to each Olson type
37683878,3879c3584
3769<     veget(:,:) = zero
3770<     frac_nobio(:,:) = zero
3771---
3772>     CASE("usda")
37733881c3586
3774<     DO vid = 1, nolson
3775---
3776>        soilclass_default=soilclass_default_usda
37773883,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"
37893893c3590
3790<     WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
3791---
3792>        ALLOCATE(textfrac_table(nscm,ntext))
37933895c3592
3794<     !   Clean up the point of the map
3795---
3796>        CALL get_soilcorr (nscm, textfrac_table)
37973899c3596
3798<        !  Let us see if all points found something in the 5km map !
3799---
3800>           ! GO through the point we have found
38013901d3597
3802<        IF ( n_found(ib) .EQ. 0 ) THEN
38033903c3599
3804<           ! Now we need to handle some exceptions
3805---
3806>           fopt = COUNT(sub_area(ib,:) > zero)
38073905,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
38153911,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
38303922a3606,3609
3831>           IF ( fopt .EQ. 0) THEN
3832>              nbexp = nbexp + 1
3833>              soilclass(ib,:) = soilclass_default
3834>              clayfraction(ib) = clayfraction_default
38353925,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))
38553943d3615
3856<           ENDIF
38573945c3617
3858<        ENDIF
3859---
3860>              !   Compute the average bare soil albedo parameters
38613946a3619
3862>              sgn = zero
38633948c3621
3864<        !  Limit the smalest vegetation fraction to 0.5%
3865---
3866>              DO ilf = 1,fopt
38673950,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
38733956,3958d3623
3874<        sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3875<        frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3876<        veget(ib,:) = veget(ib,:)/sumf
38773959a3625,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
38883963,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( cospa*coscolat + sinpa*sincolat*COS(p))
3941<        cosang(i) = cospa*coscolat + sinpa*sincolat*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>           !
39574018,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>     !
39994026c3688
4000<   !! Interpolate the Zobler soil type map
4001---
4002>   !! Calculate mean slope coef in each  model grid box from the slope map
40034028,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)
40084031,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.
40134044c3701
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
40174048,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
40234054d3708
4024<     INTEGER(i_std)               :: nbvmax
40254055a3710
4026>     REAL(r_std)  :: slope_noreinf                 ! Slope above which runoff is maximum
40274057,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
40374062d3716
4038<     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)  :: sub_area
40394064,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
40544073c3724
4055<     LOGICAL                  :: ok_interpol  ! optionnal return of aggregate_2d
4056---
4057>     LOGICAL ::           ok_interpol = .FALSE. ! optionnal return of aggregate_2d
40584078d3728
4059<     CALL get_soilcorr (classnb, textfrac_table)
40604080c3730,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
40684081a3736
4069>     slope_noreinf = 0.5
40704083,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)
40824093,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.
40934096,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)
41034105,4106d3756
4104<     ! soils_param.nc file is 1° soit texture file.
4105<     !
41064108c3758
4107<     ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
4108---
4109>     ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
41104110c3760
4111<       WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
4112---
4113>       WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
41144114c3764
4115<     ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
4116---
4117>     ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
41184116c3766
4119<       WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
4120---
4121>       WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
41224120c3770
4123<     ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4124---
4125>     ALLOCATE(slopemap(iml,jml), STAT=ALLOC_ERR)
41264122c3772
4127<       WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4128---
4129>       WRITE(numout,*) "ERROR IN ALLOCATION of slopemap : ",ALLOC_ERR
41304126c3776
4131<     ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR)
4132---
4133>     ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
41344128c3778
4135<       WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
4136---
4137>       WRITE(numout,*) "ERROR IN ALLOCATION of resol_lu : ",ALLOC_ERR
41384132,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'
41484140,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)
41564143c3789,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)
41644145c3795,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
41864149a3818,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>     !
41944153c3828
4195<           IF (soiltext(ip,jp) .GT. min_sechiba) THEN
4196---
4197>           IF (slopemap(ip,jp) .GT. min_sechiba) THEN
41984155a3831,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>           !
42214159d3855
4222<     nbvmax = 200
42234161c3857,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"
42374163,4164d3868
4238<     ok_interpol = .FALSE.
4239<     DO WHILE ( .NOT. ok_interpol )
42404169,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
42474193,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
42564202a3893,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)
42654204,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
42874223,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)))
43224255,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
43294260d3909
4330<              !
43314263,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
43444276c3913
4345<     ENDDO
4346---
4347>     WRITE(numout,*) 'Interpolation Done'
43484278,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
43554285,4288c3916
4356<     DEALLOCATE (lat_rel)
4357<     DEALLOCATE (lon_rel)
4358<     DEALLOCATE (mask)
4359<     DEALLOCATE (sub_area)
4360---
4361>     DEALLOCATE(slopemap)
43624290,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)
43724296,4298c3928
4373<   END SUBROUTINE slowproc_soilt
4374<     !
4375<
4376---
4377>   END SUBROUTINE slowproc_slope
4378_______________________________________________________________________________________________________________________
4379diff -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
4380107,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
4384129a127
4385>
4386135,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<
4482318a222,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>
4486341c247,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
4494456c366
4495<
4496---
4497>     !
4498463a374,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>         !
4611477a500,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
4621Seulement dans /login/IPSL_UTILISATEURS/MATTHIEU/VERSION_Jan_NEW/modeles/ORCHIDEE/src_sechiba: thermosoil.f90.old
4622_______________________________________________________________________________________________________________________
4623diff -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
4624369c369
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 : ', &
4628542a543,546
4629> !
4630> ! Jan : Quick fix until the steping is solved correctly in this version :-(
4631> !
4632>   for_offset = -itau_dep
4633586c590
4634<      STOP 'dim2_driver'
4635---
4636> !!     STOP 'dim2_driver'
4637886a891
4638> !
46391090a1096,1100
4640> !      petAcoef(:,:) = zero
4641> !      petBcoef(:,:) = zero
4642> !      peqAcoef(:,:) = zero
4643> !      peqBcoef(:,:) = zero
4644>       
46451197a1208
4646> !
46471657c1668,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_______________________________________________________________________________________________________________________
4653diff -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
465468c68
4655<     LOGICAL :: debug = .FALSE.
4656---
4657>     LOGICAL :: debug = .TRUE.
4658110c110,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
4663485c486
4664<   & fcontfrac
4665---
4666>   & fcontfrac, tmp_cont
4667493a495
4668>    INTEGER,DIMENSION(2) :: kloc
4669557a560,561
4670>    tmp_cont(:,:) = zero
4671> !-
4672564a569,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(u*tmp_cont)) > 100 .OR. ABS(MAXVAL(v*tmp_cont)) > 100 ) THEN
4700>       WRITE(numout,*) "Too strong wind at the surface (u or v > 100 m/s) :", &
4701>            MAXVAL(u*tmp_cont), MAXVAL(v*tmp_cont)
4702>       STOP
4703>    ENDIF
4704> !
4705678,679c714
4706< !  to bypass FPE problem on integer convertion with missing_value heigher than precision
4707<    INTEGER, PARAMETER                         :: undef_int = 999999999
4708---
4709> !
4710720c755
4711<      IF (is_watchout) THEN
4712---
4713> !!     IF (is_watchout) THEN
4714727c762
4715<      ENDIF
4716---
4717> !!     ENDIF
4718950c985
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
4722963a999
4723>                IF (check) WRITE(numout,*) "Date nm1 : ", itau_read_nm1, itau_read
47241356c1392
4725<   IF (MINVAL(tair(:,:)) < 100.) THEN
4726---
4727>   IF (MINVAL(tair(:,:)) < 100.  .AND. MINVAL(tair(:,:)) < 0.) THEN
47281381c1417
4729<            IF (tair(i,j) < 500.) THEN
4730---
4731>            IF (tair(i,j) < 500. .AND. tair(i,j) > 5.) THEN
4732_______________________________________________________________________________________________________________________
4733diff -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
47341490,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
47401671c1671
4741<       WRITE(numout,'(a1,$)') maskchar(NINT(lsm_file(i,j)))
4742---
4743>       WRITE(numout,*) maskchar(NINT(lsm_file(i,j)))
47441730c1730
4745<       WRITE(numout,'(a1,$)') maskchar(NINT(mask(i,j)))
4746---
4747>       WRITE(numout,*) maskchar(NINT(mask(i,j)))
47483058a3059
4749>   REAL,DIMENSION(1)       :: tau_tmp
47503063c3064,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, &
47553065c3067,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, &