source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_global/interpol_help.f90

Last change on this file was 4977, checked in by simon.bowring, 6 years ago

Currently running (13/02/2018) version includes all necessarily changes to include DOC in MICT code... further parametrisation necessary to equate soil pools with those of normal forcesoil restarts
  • Property svn:keywords set to HeadURL Date Author Revision
File size: 35.1 KB
Line 
1!
2! Aggregation routines. These routines allow to interpolate from the finer grid on which the
3! surface parameter is available to the coarser one of the model.
4!
5! The routines work for the fine data on a regular lat/lon grid. This grid can come in as either
6! a rank2 array or a vector. Two procedure exist which require slightly different input fields.
7!
8!
9!< $HeadURL$
10!< $Date$
11!< $Revision$
12!
13!
14MODULE interpol_help
15
16  ! Modules used :
17
18  USE constantes
19  USE mod_orchidee_para
20  USE grid
21  USE interregxy
22  USE mod_orchidee_transfert_para
23
24  IMPLICIT NONE
25
26  PRIVATE
27  PUBLIC aggregate, aggregate_p
28  !
29  INTERFACE aggregate
30     MODULE PROCEDURE aggregate_2d, aggregate_vec
31  END INTERFACE
32  !
33  INTERFACE aggregate_p
34     MODULE PROCEDURE aggregate_2d_p, aggregate_vec_p
35  END INTERFACE
36  !
37  LOGICAL, PARAMETER                              :: check_grid=.FALSE.
38  !
39CONTAINS
40  !
41  ! This routing will get for each point of the coarse grid the
42  ! indexes of the finer grid and the area of overlap.
43  ! This routine is designed for a fine grid which is regular in lat/lon.
44  !
45  SUBROUTINE aggregate_2d (nbpt, lalo, neighbours, resolution, contfrac, &
46       &                iml, jml, lon_rel, lat_rel, mask, callsign, &
47       &                incmax, indinc, areaoverlap, ok, opt_nbpt_start, opt_nbpt_end)
48
49    USE grid, ONLY : global
50
51    !
52    ! INPUT
53    !
54    INTEGER(i_std), INTENT(in)   :: nbpt                 ! Number of points for which the data needs to be interpolated
55    REAL(r_std), INTENT(in)      :: lalo(nbpt,2)         ! Vector of latitude and longitudes (beware of the order !)
56    INTEGER(i_std), INTENT(in)   :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
57    REAL(r_std), INTENT(in)      :: resolution(nbpt,2)   ! The size in km of each grid-box in X and Y
58    REAL(r_std), INTENT(in)      :: contfrac(nbpt)       ! Fraction of land in each grid box.
59    INTEGER(i_std), INTENT(in)   :: iml, jml             ! Size of the finer grid
60    REAL(r_std), INTENT(in)      :: lon_rel(iml, jml)    ! Longitudes for the finer grid
61    REAL(r_std), INTENT(in)      :: lat_rel(iml, jml)    ! Latitudes for the finer grid
62    INTEGER(i_std), INTENT(in)   :: mask(iml, jml)       ! Mask which retains only the significative points
63                                                         ! of the fine grid.
64    CHARACTER(LEN=*), INTENT(in) :: callsign             ! Allows to specify which variable is beeing treated
65    INTEGER(i_std), INTENT(in)   :: incmax              ! Maximum point of the fine grid we can store.
66    INTEGER(i_std), OPTIONAL, INTENT(in)    :: opt_nbpt_start            ! Input Start grid cell interpolation
67    INTEGER(i_std), OPTIONAL, INTENT(in)    :: opt_nbpt_end              ! Input End grid cell interpolation
68    !
69    ! Output
70    !
71    INTEGER(i_std), INTENT(out)  :: indinc(:,:,:)
72    REAL(r_std), INTENT(out)      :: areaoverlap(:,:)
73    LOGICAL, OPTIONAL, INTENT(out)      :: ok            ! return code
74    !
75    ! Local Variables
76    !
77    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
78    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel
79    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: searchind
80    REAL(r_std) :: lon_up, lon_low, lat_up, lat_low, lat_resol, lon_resol
81    REAL(r_std) :: coslat, ax, ay, sgn, lonrel, lolowrel, louprel
82    INTEGER(i_std) :: fopt, fopt_max, ip, jp, ib, i, itmp, iprog, nbind
83    REAL(r_std) :: domain_minlon,domain_maxlon,domain_minlat,domain_maxlat
84    INTEGER(i_std) :: minLon(1), maxLon(1)
85    INTEGER(i_std) :: nbpt_start            ! Start grid cell interpolation
86    INTEGER(i_std) :: nbpt_end              ! End grid cell interpolation
87    INTEGER(i_std) :: landpoint_idx 
88
89    INTEGER                  :: ALLOC_ERR
90    LOGICAL :: err_fopt
91    err_fopt = .FALSE.
92    !
93    ! Some inital assignmens
94    !
95    areaoverlap(:,:) = moins_un
96    indinc(:,:,:) = zero
97
98    ALLOCATE (laup_rel(iml,jml), STAT=ALLOC_ERR)
99    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of laup_rel','','')
100
101    ALLOCATE (loup_rel(iml,jml), STAT=ALLOC_ERR)
102    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of loup_rel','','')
103
104    ALLOCATE (lalow_rel(iml,jml), STAT=ALLOC_ERR)
105    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lalow_rel','','')
106
107    ALLOCATE (lolow_rel(iml,jml), STAT=ALLOC_ERR)
108    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lolow_rel','','')
109
110    ALLOCATE (lat_ful(iml+2,jml+2), STAT=ALLOC_ERR)
111    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lat_ful','','')
112
113    ALLOCATE (lon_ful(iml+2,jml+2), STAT=ALLOC_ERR)
114    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lon_ful','','')
115
116    ALLOCATE (searchind(iml*jml,2), STAT=ALLOC_ERR)
117    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of searchind','','')
118
119    nbpt_start = 1
120    nbpt_end = nbpt
121    IF (PRESENT(opt_nbpt_start) .AND. PRESENT(opt_nbpt_end)) THEN
122        nbpt_start = opt_nbpt_start
123        nbpt_end = opt_nbpt_end
124    ENDIF
125
126    IF (PRESENT(ok)) ok = .TRUE.
127    !
128    !    Duplicate the border assuming we have a global grid going from west to east
129    !
130    lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml)
131    lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml)
132    !
133    IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN
134       lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)
135       lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
136    ELSE
137       lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360
138       lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
139    ENDIF
140
141    IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN
142       lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)
143       lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
144    ELSE
145       lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360
146       lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
147    ENDIF
148    !
149    lat_resol = lat_rel(1,2) - lat_rel(1,1)
150    lat_ful(2:iml+1,1) = lat_rel(1,1) - lat_resol
151    lat_ful(2:iml+1,jml+2) = lat_rel(1,jml) + lat_resol
152    lat_ful(1,1) = lat_ful(iml+1,1)
153    lat_ful(iml+2,1) = lat_ful(2,1)
154    lat_ful(1,jml+2) = lat_ful(iml+1,jml+2)
155    lat_ful(iml+2,jml+2) = lat_ful(2,jml+2)
156    !
157    ! Add the longitude lines to the top and bottom
158    !
159    lon_resol = lon_ful(3,2) - lon_ful(2,2)
160    lon_ful(:,1) = lon_ful(:,2) 
161    lon_ful(:,jml+2) = lon_ful(:,jml+1) 
162    lon_ful(1,:) = lon_ful(2,1) - lon_resol
163    lon_ful(iml+2,:) = lon_ful(iml+1,1) + lon_resol
164    !
165    !  Get the upper and lower limits of each grid box
166    !
167    DO ip=1,iml
168       DO jp=1,jml
169          !
170          loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),&
171               & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
172          lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),&
173               & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
174          laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),&
175               & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
176          lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),&
177               & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
178          !
179       ENDDO
180    ENDDO
181    IF (check_grid) THEN
182       WRITE(numout,*) "================================"
183       WRITE(numout,*) "interpol_aggregate_2d : "
184       WRITE(numout,*) "lalo(:,1) :",lalo(:,1)
185       WRITE(numout,*) "lalo(:,2) :",lalo(:,2)
186       WRITE(numout,*) "Map meshes : "
187       WRITE(numout,*) "lat read(1,:) :",lat_rel(1,:)
188       WRITE(numout,*) "lat_ful(1,:) :",lat_ful(1,:)
189       WRITE(numout,*) "lat_ful(2,:) :",lat_ful(2,:)
190       WRITE(numout,*) "lalow_rel(1,:) :",lalow_rel(1,:)
191       WRITE(numout,*) "laup_rel(1,:) :",laup_rel(1,:)
192       WRITE(numout,*) "================================"
193       WRITE(numout,*) "lon read(:,1) :",lon_rel(:,1)
194       WRITE(numout,*) "lon_ful(:,1) :",lon_ful(:,1)
195       WRITE(numout,*) "lon_ful(:,2) :",lon_ful(:,2)
196       WRITE(numout,*) "lolow_rel(:,1) :",lolow_rel(:,1)
197       WRITE(numout,*) "loup_rel(:,1) :",loup_rel(:,1)
198       WRITE(numout,*) "================================"
199    ENDIF
200    !
201    !
202    !  To speedup calculations we will get the limits of the domain of the
203    !  coarse grid and select all the points of the fine grid which are potentialy
204    !  in this domain.
205    !
206    !
207    minLon = MINLOC(lalo(1:nbpt,2))
208    coslat = MAX(COS(lalo(minLon(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
209    domain_minlon = lalo(minLon(1),2) - resolution(minLon(1),1)/(2.0*coslat)
210    maxLon = MAXLOC(lalo(1:nbpt,2))
211    coslat = MAX(COS(lalo(maxLon(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
212    domain_maxlon = lalo(maxLon(1),2) + resolution(maxLon(1),1)/(2.0*coslat)
213    !
214    coslat = pi/180. * R_Earth
215    domain_minlat = MINVAL(lalo(1:nbpt,1)) - resolution(maxLon(1),2)/(2.0*coslat)
216    domain_maxlat = MAXVAL(lalo(1:nbpt,1)) + resolution(maxLon(1),2)/(2.0*coslat)
217    !
218    IF (check_grid) THEN
219       WRITE(numout,*) "indices min/max of longitude :",minLon,maxLon, &
220            & "; longitude min/max : ",lalo(minLon,1),lalo(maxLon,1)
221       WRITE(numout,*) "Domain for coarse grid :"
222       WRITE(numout,*) '(',domain_minlat,',',domain_minlon,')',&
223   &                   '(',domain_maxlat,',',domain_maxlon,')'
224       WRITE(numout,*) "================================"
225    ENDIF
226    !
227    ! we list a first approximation of all point we will need to
228    ! scan to fill our coarse grid.
229    !
230    IF ( global ) THEN
231       ! Here we do the entire globe
232       WRITE(numout,*) 'In aggregate_p : do interpolation to global model domain'
233       nbind=0
234       DO jp=1,jml
235          DO ip=1,iml
236             IF (mask(ip,jp) == 1 ) THEN
237                nbind = nbind + 1
238                searchind(nbind,1) = ip
239                searchind(nbind,2) = jp
240             ENDIF
241          ENDDO
242       ENDDO
243       !
244    ELSE
245       ! Now we get a limited number of points
246       WRITE(numout,*) 'In aggregate_p : do interpolation to regional model domain'
247       nbind=0
248       DO jp=1,jml
249          DO ip=1,iml
250             IF ( loup_rel(ip,jp) >= domain_minlon .AND. lolow_rel(ip,jp) <= domain_maxlon .AND.&
251               &  laup_rel(ip,jp) >= domain_minlat .AND. lalow_rel(ip,jp) <= domain_maxlat ) THEN
252                IF (mask(ip,jp) == 1 ) THEN
253                   nbind = nbind + 1
254                   searchind(nbind,1) = ip
255                   searchind(nbind,2) = jp
256                ENDIF
257             ENDIF
258          ENDDO
259       ENDDO
260    ENDIF
261    !
262    WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid'
263    !
264    WRITE(numout,*) 'Aggregate_2d : ', callsign
265#ifdef INTERPOL_ADVANCE
266    WRITE(numout,'(2a40)')'0%--------------------------------------', &
267         & '------------------------------------100%'
268#endif
269    !
270    !   Now we take each grid point and find out which values from the forcing we need to average
271    !
272    fopt_max = -1
273    DO ib = nbpt_start, nbpt_end
274       landpoint_idx = ib - nbpt_start + 1
275       !
276       !   Give a progress meter
277       !
278#ifdef INTERPOL_ADVANCE
279       iprog = NINT(REAL(ib,r_std)/REAL(nbpt,r_std)*79.) - NINT(REAL(ib-1,r_std)/REAL(nbpt,r_std)*79.)
280       IF ( iprog .NE. 0 ) THEN
281          WRITE(numout,'(a1,$)') 'x'
282       ENDIF
283#endif
284       !
285       !  We find the 4 limits of the grid-box. As we transform the resolution of the model
286       !  into longitudes and latitudes we do not have the problem of periodicity.
287       !  coslat is a help variable here !
288       !
289       coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth
290       !
291       lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) 
292       lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) 
293       !
294       coslat = pi/180. * R_Earth
295       !
296       lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) 
297       lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) 
298       !
299       !  Find the grid boxes from the data that go into the model's boxes
300       !  We still work as if we had a regular grid ! Well it needs to be localy regular so
301       !  so that the longitude at the latitude of the last found point is close to the one
302       !  of the next point.
303       !
304       fopt = zero
305       !
306       DO i=1,nbind
307          !
308          ip = searchind(i,1)
309          jp = searchind(i,2)
310          !
311          !  Either the center of the data grid point is in the interval of the model grid or
312          !  the East and West limits of the data grid point are on either sides of the border of
313          !  the data grid.
314          !
315          !  To do that correctly we have to check if the grid box sits on the date-line.
316          !
317          IF (lon_low >= -180 .AND. lon_up <= 180.0) THEN
318             lonrel = lon_rel(ip,jp)
319             lolowrel = lolow_rel(ip,jp)
320             louprel = loup_rel(ip,jp)
321          ELSE IF ( lon_low < -180.0 ) THEN
322             ! -179 -> -179
323             ! 179 -> -181
324             lonrel = MOD( lon_rel(ip,jp) - 360.0, 360.0)
325             lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
326             louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
327             !
328          ELSE IF ( lon_up > 180.0 ) THEN
329             ! -179 -> 181
330             !  179 -> 179
331             lonrel = MOD( lon_rel(ip,jp) + 360., 360.0)
332             lolowrel = MOD( lolow_rel(ip,jp) + 360., 360.0)
333             louprel = MOD( loup_rel(ip,jp) + 360., 360.0)
334          ENDIF
335          !
336          !
337          !
338          IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
339               & lolowrel < lon_low .AND.  louprel > lon_low .OR. &
340               & lolowrel < lon_up  .AND.  louprel > lon_up ) THEN
341             !
342             ! Now that we have the longitude let us find the latitude
343             !
344             IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. &
345                  & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.&
346                  & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN
347                   !
348                fopt = fopt + 1
349                IF ( fopt > incmax) THEN
350                   err_fopt=.TRUE.
351                   EXIT
352                ELSE
353                   !
354                   ! If we sit on the date line we need to do the same transformations as above.
355                   !
356                   IF ( lon_low >= -180.0 .AND. lon_up <= 180.0 ) THEN
357                      lolowrel = lolow_rel(ip,jp)
358                      louprel = loup_rel(ip,jp)
359                   ELSE IF ( lon_low < -180.0 ) THEN
360                      lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
361                      louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
362                      !
363                   ELSE IF ( lon_up > 180.0 ) THEN
364                      lolowrel = MOD( lolow_rel(ip,jp) + 360., 360.0)
365                      louprel = MOD( loup_rel(ip,jp) + 360., 360.0)
366!                   ELSE
367!                      lolowrel = lolow_rel(ip,jp)
368!                      louprel = loup_rel(ip,jp)
369                   ENDIF
370                   !
371                   ! Get the area of the fine grid in the model grid
372                   !
373                   coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )
374                   ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat
375                   ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth
376                   !
377                   areaoverlap(landpoint_idx, fopt) = ax*ay
378                   indinc(landpoint_idx, fopt, 1) = ip
379                   indinc(landpoint_idx, fopt, 2) = jp
380                   !
381                   ! If this point was 100% within the grid then we can de-select it from our
382                   ! list as it can not be in another mesh of the coarse grid.
383                   !
384                   IF ( louprel < lon_up .AND. lolowrel > lon_low .AND. &
385                     &  laup_rel(ip,jp) < lat_up .AND. lalow_rel(ip,jp) > lat_low ) THEN
386                      searchind(i,1) = 0
387                      searchind(i,2) = 0
388                   ENDIF
389                   !
390                ENDIF
391             ENDIF       ! IF lat
392          ENDIF          ! IF lon
393       ENDDO
394
395       IF (err_fopt) THEN
396          WRITE(numout,*) 'Working on variable :', callsign
397          WRITE(numout,*) 'Reached value ', fopt,' for fopt on point', ib, lalo(ib,2), lalo(ib,1)
398          CALL ipslerr_p(2,'aggregate_2d', &
399               'Working on variable :'//callsign, &
400               'Reached incmax value for fopt.',&
401               'Please increase incmax in subroutine calling aggregate')                   
402          IF (PRESENT(ok)) THEN
403             ok = .FALSE.
404             RETURN
405          ELSE
406             CALL ipslerr_p(3,'aggregate_2d','Stop now','','')
407          ENDIF
408       ENDIF
409       fopt_max = MAX ( fopt, fopt_max )
410       !
411       ! De-select the marked points
412       !
413       itmp = nbind
414       nbind = 0
415       DO i=1,itmp
416          IF ( searchind(i,1) > 0 .AND. searchind(i,2) > 0 ) THEN
417             nbind = nbind + 1
418             searchind(nbind,1) = searchind(i,1)
419             searchind(nbind,2) = searchind(i,2)
420          ENDIF
421       ENDDO
422       !
423    ENDDO
424    !
425    DO ib=nbpt_start, nbpt_end
426       DO fopt=1,incmax
427          IF (( indinc(landpoint_idx,fopt,1) == 0 .AND. indinc(landpoint_idx,fopt,2) > 0) .OR.&
428               & ( indinc(landpoint_idx,fopt,2) == 0 .AND. indinc(landpoint_idx,fopt,1) > 0) ) THEN
429             WRITE(*,*) "aggregate_2d PROBLEM : point =",ib, fopt," Indicies = ", &
430                  & indinc(landpoint_idx,fopt,1), indinc(landpoint_idx,fopt,2), areaoverlap(landpoint_idx,fopt)
431          ENDIF
432       ENDDO
433    ENDDO
434
435
436    WRITE(numout,*) ""
437    WRITE(numout,*) "aggregate_2D nbvmax = ",incmax, "max used = ",fopt_max
438    !
439    ! Do some memory management.
440    !
441    DEALLOCATE (laup_rel)
442    DEALLOCATE (loup_rel)
443    DEALLOCATE (lalow_rel)
444    DEALLOCATE (lolow_rel)
445    DEALLOCATE (lat_ful)
446    DEALLOCATE (lon_ful)
447    DEALLOCATE (searchind)
448    !
449    ! Close the progress meter
450    !
451    WRITE(numout,*) '    '
452    !
453  END SUBROUTINE aggregate_2d
454
455  !
456  ! This routing will get for each point of the coarse grid the
457  ! indexes of the finer grid and the area of overlap.
458  ! This routine is designed for a fine grid which is regular in meters along lat lon axes.
459  !
460  SUBROUTINE aggregate_vec (nbpt, lalo, neighbours, resolution, contfrac, &
461       &                iml, lon_rel, lat_rel, resol_lon, resol_lat, callsign, &
462       &                incmax, indinc, areaoverlap, ok)
463    !
464    ! INPUT
465    !
466    INTEGER(i_std), INTENT(in)   :: nbpt                 ! Number of points for which the data needs to be interpolated
467    REAL(r_std), INTENT(in)      :: lalo(nbpt,2)         ! Vector of latitude and longitudes (beware of the order !)
468    INTEGER(i_std), INTENT(in)   :: neighbours(nbpt,NbNeighb)   ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
469    REAL(r_std), INTENT(in)      :: resolution(nbpt,2)   ! The size in km of each grid-box in X and Y
470    REAL(r_std), INTENT(in)      :: contfrac(nbpt)       ! Fraction of land in each grid box.
471    INTEGER(i_std), INTENT(in)   :: iml                  ! Size of the finer grid
472    REAL(r_std), INTENT(in)      :: lon_rel(iml)         ! Longitudes for the finer grid
473    REAL(r_std), INTENT(in)      :: lat_rel(iml)         ! Latitudes for the finer grid
474    REAL(r_std), INTENT(in)      :: resol_lon, resol_lat ! Resolution in meters of the fine grid
475    CHARACTER(LEN=*), INTENT(in) :: callsign             ! Allows to specify which variable is beeing treated
476    INTEGER(i_std), INTENT(in)   :: incmax              ! Maximum point of the fine grid we can store.
477    !
478    ! Output
479    !
480    INTEGER(i_std), INTENT(out)  :: indinc(nbpt,incmax)
481    REAL(r_std), INTENT(out)      :: areaoverlap(nbpt,incmax)
482    LOGICAL, OPTIONAL, INTENT(out)      :: ok            ! return code
483    !
484    ! Local Variables
485    !
486    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: searchind
487    REAL(r_std) :: lon_up, lon_low, lat_up, lat_low
488    REAL(r_std) :: coslat, ax, ay, lonrel, lolowrel, louprel
489    REAL(r_std) :: latrel, lauprel, lalowrel
490    INTEGER(i_std), DIMENSION(nbpt) :: fopt
491    INTEGER(i_std) :: fopt_max, not_found_fopt
492    INTEGER(i_std) :: ip, ib, i, j, itmp, iprog, nbind, pp, ipp
493    REAL(r_std) :: domain_minlon,domain_maxlon,domain_minlat,domain_maxlat
494    REAL(r_std) :: minlon, minlat, mini
495    INTEGER(i_std) :: ff(1), incp
496    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: fine_ind
497    INTEGER(i_std) :: pos_pnt(5)
498    INTEGER                  :: ALLOC_ERR
499    !
500    LOGICAL :: err_fopt
501    err_fopt = .FALSE.
502    !
503    ! Some inital assignmens
504    !
505    areaoverlap(:,:) = moins_un
506    indinc(:,:) = zero
507
508    ALLOCATE (searchind(iml), STAT=ALLOC_ERR)
509    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_vec','ERROR IN ALLOCATION of searchbind','','')
510
511    IF (PRESENT(ok)) ok = .TRUE.
512    !
513    !  To speedup calculations we will get the limits of the domain of the
514    !  coarse grid and select all the points of the fine grid which are potentialy
515    !  in this domain.
516    !
517    !
518    ff = MINLOC(lalo(:,2))
519    coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
520    domain_minlon = lalo(ff(1),2) - resolution(ff(1),1)/(2.0*coslat)
521    ff = MAXLOC(lalo(:,2))
522    coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
523    domain_maxlon = lalo(ff(1),2) + resolution(ff(1),1)/(2.0*coslat)
524    !
525    coslat = pi/180. * R_Earth
526    domain_minlat = MINVAL(lalo(:,1)) - resolution(ff(1),2)/(2.0*coslat)
527    domain_maxlat = MAXVAL(lalo(:,1)) + resolution(ff(1),2)/(2.0*coslat)
528    !
529    ! Find appropriate resolution for index table
530    !
531    ff=MINLOC(resolution(:,1))
532    coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
533    minlon=resolution(ff(1),1)/(2.0*coslat)
534    ff=MINLOC(resolution(:,2))
535    coslat = pi/180. * R_Earth
536    minlat=resolution(ff(1),2)/(2.0*coslat)
537    mini=MIN(minlon, minlat)
538    !
539    ! This interpolation only works if the model grid is coarser than the data grid
540    !
541    IF (MINVAL(resolution(:,1)) < resol_lon .OR. MINVAL(resolution(:,2)) < resol_lat) THEN
542       WRITE(numout,*) " === WARNING == "
543       WRITE(numout,*) "Resolution minima of the model (lon, lat) : ", &
544            & MINVAL(resolution(:,1)), MINVAL(resolution(:,2))
545       WRITE(numout,*) "Resolution of the file to be interpolated (fine grid) : ", resol_lon, resol_lat
546       WRITE(numout,*) "This interpolation assumes that we aggregate from a fine grid to a coarser grid"
547       WRITE(numout,*) "In the data submitted it apears that the model is runing on a finer grid than the data"
548    ENDIF
549    !
550    incp = 10
551    IF (mini < 0.1) THEN
552       incp=100
553    ELSE IF (mini < 0.01) THEN
554       incp = 1000
555    ENDIF
556    !
557    ! Allocate the needed memory for fine_ind
558    !
559    ALLOCATE (fine_ind(NINT(domain_minlon*incp)-2:NINT(domain_maxlon*incp)+2, &
560         & NINT(domain_minlat*incp)-2:NINT(domain_maxlat*incp)+2), STAT=ALLOC_ERR)
561    IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_vec','ERROR IN ALLOCATION of find_ind','','')
562    !
563    ! Generate a quick access table for the coarse grid
564    !
565    fine_ind(:,:) = zero
566    !
567    DO ib=1,nbpt
568       coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth
569       !
570       lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) 
571       lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) 
572       !
573       coslat = pi/180. * R_Earth
574       !
575       lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) 
576       lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) 
577       !
578       fine_ind(NINT(lon_low*incp):NINT(lon_up*incp),NINT(lat_low*incp):NINT(lat_up*incp))=ib
579       !
580    ENDDO
581    !
582    WRITE(numout,*) 'Domaine LON range : ', domain_minlon, domain_maxlon
583    WRITE(numout,*) 'Domaine LAT range : ', domain_minlat, domain_maxlat
584    !
585    ! we list a first approximation of all point we will need to
586    ! scan to fill our coarse grid.
587    !
588    IF ( domain_minlon <= -179.5 .AND. domain_maxlon >= 179.5 .AND. &
589      &  domain_minlat <= -89.5  .AND. domain_maxlat >= 89.5 ) THEN
590       ! Here we do the entire globe
591       nbind=0
592       DO ip=1,iml
593          nbind = nbind + 1
594          searchind(nbind) = ip
595       ENDDO
596       !
597    ELSE
598       ! Now we get a limited number of points
599       nbind=0
600       DO ip=1,iml
601          ! Compute the limits of the meshes of the fine grid
602          coslat = MAX(COS(lat_rel(ip) * pi/180. ), mincos )*pi/180. * R_Earth
603          louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), 180.)
604          lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), -180.)
605          coslat = pi/180. * R_Earth
606          lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), 90.)
607          lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), -90.)
608          !
609          IF ( louprel >= domain_minlon .AND. lolowrel <= domain_maxlon .AND.&
610            &  lauprel >= domain_minlat .AND. lalowrel <= domain_maxlat ) THEN
611             nbind = nbind + 1
612             searchind(nbind) = ip
613          ENDIF
614       ENDDO
615    ENDIF
616    !
617    WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid and ', nbpt, 'for the coarse grid'
618    !
619    WRITE(numout,*) 'Aggregate_vec : ', callsign
620    !
621    !   Now we take each grid point and find out which values from the forcing we need to average
622    !
623    fopt(:) = zero
624    fopt_max = -1
625    !
626    !
627    !
628    loopnbind : DO i=1,nbind
629       !
630       !
631       ip = searchind(i)
632       !
633       !  Either the center of the data grid point is in the interval of the model grid or
634       !  the East and West limits of the data grid point are on either sides of the border of
635       !  the data grid.
636       !
637       lonrel = lon_rel(ip)
638       coslat = MAX(COS(lat_rel(ip) * pi/180. ), mincos )*pi/180. * R_Earth
639       louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), domain_maxlon)
640       lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), domain_minlon)
641       !
642       latrel = lat_rel(ip)
643       coslat = pi/180. * R_Earth
644       lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), domain_maxlat)
645       lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), domain_minlat)
646       !
647       !
648       pos_pnt(:) = zero
649       ipp = zero
650       pp = fine_ind(NINT(lonrel*incp),NINT(latrel*incp))
651       !
652       IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
653          pos_pnt(ipp+1) = pp
654          ipp = ipp + 1
655       ENDIF
656       pp = fine_ind(NINT(louprel*incp),NINT(lauprel*incp))
657       !
658       IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
659          pos_pnt(ipp+1) = pp
660          ipp = ipp + 1
661       ENDIF
662       pp = fine_ind(NINT(louprel*incp),NINT(lalowrel*incp))
663       !
664       IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
665          pos_pnt(ipp+1) = pp
666          ipp = ipp + 1
667       ENDIF
668       pp = fine_ind(NINT(lolowrel*incp),NINT(lauprel*incp))
669       !
670       IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
671          pos_pnt(ipp+1) = pp
672          ipp = ipp + 1
673       ENDIF
674       pp = fine_ind(NINT(lolowrel*incp),NINT(lalowrel*incp))
675       !
676       IF (COUNT(pos_pnt(:) == pp) == zero ) THEN
677          pos_pnt(ipp+1) = pp
678          ipp = ipp + 1
679       ENDIF
680       !
681       !
682       IF ( ipp > zero ) THEN
683          !
684          DO pp=1,ipp
685             ib = pos_pnt(pp)
686             !
687             !  We find the 4 limits of the grid-box. As we transform the resolution of the model
688             !  into longitudes and latitudes we do not have the problem of periodicity.
689             !  coslat is a help variable here !
690             !
691             coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth
692             !
693             lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) 
694             lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) 
695             !
696             coslat = pi/180. * R_Earth
697             !
698             lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) 
699             lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) 
700             !
701             IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
702                  & lolowrel < lon_low .AND.  louprel > lon_low .OR. &
703                  & lolowrel < lon_up  .AND.  louprel > lon_up ) THEN
704                !
705                ! Now that we have the longitude let us find the latitude
706                !             
707                IF ( latrel > lat_low .AND. latrel < lat_up .OR. &
708                     & lalowrel < lat_low .AND. lauprel > lat_low .OR.&
709                     & lalowrel < lat_up .AND. lauprel > lat_up) THEN
710                   !
711                   fopt(ib) = fopt(ib) + 1
712                   fopt_max = MAX ( fopt(ib), fopt_max )
713                   !
714                   IF ( fopt(ib) > incmax) THEN
715                      err_fopt=.TRUE.
716                      EXIT loopnbind
717                   ELSE
718                      !
719                      ! Get the area of the fine grid in the model grid
720                      !
721                      coslat = MAX( COS( lat_rel(ip) * pi/180. ), mincos )
722                      ax = (MIN(lon_up,louprel)-MAX(lon_low,lolowrel))*pi/180. * R_Earth * coslat
723                      ay = (MIN(lat_up,lauprel)-MAX(lat_low,lalowrel))*pi/180. * R_Earth
724                      !
725                      areaoverlap(ib, fopt(ib)) = ax*ay
726                      indinc(ib, fopt(ib)) = ip
727                      !
728                   ENDIF
729                ENDIF
730             ENDIF
731          ENDDO
732       ENDIF
733    ENDDO loopnbind
734    !
735    IF (err_fopt) THEN
736       WRITE(numout,*) 'Reached value ', fopt(ib),' for fopt on point', ib
737       CALL ipslerr_p(2,'aggregate_vec (nbpt < nbind)', &
738            'Working on variable :'//callsign, &
739            'Reached incmax value for fopt.',&
740            'Please increase incmax in subroutine calling aggregate')
741       IF (PRESENT(ok)) THEN
742          ok = .FALSE.
743          RETURN
744       ELSE
745          CALL ipslerr_p(3,'aggregate_vec','Stop now','','')
746       ENDIF
747    ENDIF
748    !
749    WRITE(numout,*) 
750    not_found_fopt = COUNT(fopt(:) .EQ. zero)
751    WRITE(numout,*) "aggregate_vec : ",not_found_fopt, &
752         & "did not find any corresponding data in the input file."
753    WRITE(numout,*) "aggregate_vec : This is ", not_found_fopt/REAL(nbpt, r_std)*100., &
754         & " % of the grid"
755    WRITE(numout,*) "aggregate_vec : nbvmax = ",incmax, "max used = ",fopt_max
756    !
757    ! Do some memory management.
758    !
759    DEALLOCATE (searchind)
760    DEALLOCATE (fine_ind)
761    !
762    ! Close the progress meter
763    !
764    WRITE(numout,*) '    '
765    !
766  END SUBROUTINE aggregate_vec
767!
768!
769  SUBROUTINE aggregate_vec_p(nbpt, lalo, neighbours, resolution, contfrac,          &
770       &                 iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
771       &                 nbvmax, sub_index, sub_area, ok)
772   
773    IMPLICIT NONE
774   
775    INTEGER(i_std), INTENT(in)   :: nbpt                 
776    REAL(r_std), INTENT(in)      :: lalo(nbpt,2)
777    INTEGER(i_std), INTENT(in)   :: neighbours(nbpt,NbNeighb)   
778    REAL(r_std), INTENT(in)      :: resolution(nbpt,2)   
779    REAL(r_std), INTENT(in)      :: contfrac(nbpt)       
780    INTEGER(i_std), INTENT(in)   :: iml                 
781    REAL(r_std), INTENT(in)      :: lon_ful(iml)         
782    REAL(r_std), INTENT(in)      :: lat_ful(iml)         
783    REAL(r_std), INTENT(in)      :: resol_lon, resol_lat 
784    CHARACTER(LEN=*), INTENT(in) :: callsign             
785    INTEGER(i_std), INTENT(in)   :: nbvmax             
786    INTEGER(i_std), INTENT(out)  :: sub_index(nbpt,nbvmax)
787    REAL(r_std), INTENT(out)     :: sub_area(nbpt,nbvmax) 
788    LOGICAL, OPTIONAL, INTENT(out)  :: ok            ! return code
789
790    INTEGER(i_std)  :: sub_index_g(nbp_glo,nbvmax)
791    REAL(r_std)       :: sub_area_g(nbp_glo,nbvmax)
792
793    IF ( INDEX(GridType,"RegLonLat") > 0) THEN
794       IF (is_root_prc) CALL aggregate_vec(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
795            &                                  iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign,   &
796            &                                  nbvmax, sub_index_g, sub_area_g, ok)
797    ELSE IF ( INDEX(GridType,"RegXY") > 0) THEN
798       IF ( proj_stack(1)%code > undef_int-1 ) THEN
799          CALL ipslerr(3, "aggregate_vec_p", "RegXY projection was not intialized.", &
800               &       "proj_stack(1)%code is undefined","")
801       ENDIF
802       IF (is_root_prc) CALL interregxy_aggrve(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
803            &                                  iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign,   &
804            &                                  nbvmax, sub_index_g, sub_area_g, ok)
805    ELSE
806       CALL ipslerr(3, "aggregate_vec_p", "Interpolation is only possible for regular lat/lon grids for the moment.", &
807            &       "More general interpolation will be implemented soon.","")
808    ENDIF
809    !
810    CALL BCAST(ok)
811    CALL scatter(sub_index_g,sub_index)
812    CALL scatter(sub_area_g,sub_area)
813   
814   
815  END SUBROUTINE aggregate_vec_p
816
817  SUBROUTINE aggregate_2d_p(nbpt, lalo, neighbours, resolution, contfrac,          &
818       &                 iml, jml, lon_ful, lat_ful, mask, callsign, &
819       &                 nbvmax, sub_index, sub_area, ok)
820   
821    IMPLICIT NONE
822   
823    INTEGER(i_std), INTENT(in)     :: nbpt                 
824    REAL(r_std), INTENT(in)        :: lalo(nbpt,2)
825    INTEGER(i_std), INTENT(in)     :: neighbours(nbpt,NbNeighb)   
826    REAL(r_std), INTENT(in)        :: resolution(nbpt,2)   
827    REAL(r_std), INTENT(in)        :: contfrac(nbpt)       
828    INTEGER(i_std), INTENT(in)     :: iml,jml                 
829    REAL(r_std), INTENT(in)        :: lon_ful(iml,jml)         
830    REAL(r_std), INTENT(in)        :: lat_ful(iml,jml)         
831    INTEGER(i_std), INTENT(in)     :: mask(iml, jml)
832    CHARACTER(LEN=*), INTENT(in)   :: callsign             
833    INTEGER(i_std), INTENT(in)     :: nbvmax             
834    INTEGER(i_std), INTENT(out)    :: sub_index(nbpt,nbvmax,2)
835    REAL(r_std), INTENT(out)       :: sub_area(nbpt,nbvmax) 
836    LOGICAL, OPTIONAL, INTENT(out) :: ok            ! return code
837    INTEGER(i_std)                 :: nbp_start, nbp_end
838
839
840    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index_g 
841    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: sub_area_g 
842
843    IF ( INDEX(GridType,"RegLonLat") > 0) THEN
844       nbp_start = nbp_mpi_para_begin(mpi_rank)
845       nbp_end = nbp_mpi_para_end(mpi_rank)
846
847        CALL aggregate_2d(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
848       &                   iml, jml, lon_ful, lat_ful, mask, callsign,   &
849       &                   nbvmax, sub_index, sub_area, ok,          &
850       &                   nbp_start, nbp_end)
851    ELSE IF ( INDEX(GridType,"RegXY") > 0) THEN
852        ALLOCATE(sub_index_g(nbp_glo, nbvmax, 2))
853        ALLOCATE(sub_area_g(nbp_glo, nbvmax))
854
855       IF ( proj_stack(1)%code > undef_int-1 ) THEN
856          CALL ipslerr(3, "aggregate_2d_p", "RegXY projection was not intialized.", &
857               &       "proj_stack(1)%code is undefined","")
858       ENDIF
859       IF (is_root_prc) CALL interregxy_aggr2d(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
860            &                               iml, jml, lon_ful, lat_ful, mask, callsign,   &
861            &                               nbvmax, sub_index_g, sub_area_g, ok)
862
863       CALL BCAST(ok)
864       CALL scatter(sub_index_g,sub_index)
865       CALL scatter(sub_area_g,sub_area)
866
867       DEALLOCATE(sub_index_g)
868       DEALLOCATE(sub_area_g)
869    ELSE
870       CALL ipslerr(3, "aggregate_2d_p", "Interpolation is only possible for regular lat/lon grids for the moment.", &
871            &       "More general interpolation will be implemented soon.","")
872
873    ENDIF
874    !
875   
876  END SUBROUTINE aggregate_2d_p
877!
878END MODULE interpol_help
Note: See TracBrowser for help on using the repository browser.