source: branches/publications/ORCHIDEE-MICT-BIOENERGY_r7298/src_global/interpol_help.f90 @ 7299

Last change on this file since 7299 was 7297, checked in by wei.li, 3 years ago

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