source: branches/publications/ORCHIDEE-PEAT_r5488/src_global/interpol_help.f90 @ 7746

Last change on this file since 7746 was 4806, checked in by chunjing.qiu, 7 years ago

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