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

Last change on this file since 7746 was 4693, checked in by josefine.ghattas, 7 years ago

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