source: branches/publications/ORCHIDEE_gmd-2018-57/src_global/interpol_help.f90 @ 5143

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