source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/interpol_help.f90 @ 7683

Last change on this file since 7683 was 3447, checked in by josefine.ghattas, 8 years ago

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