source: branches/publications/ORCHIDEE-MUSLE-r6129/src_global/interpol_help.f90 @ 7346

Last change on this file since 7346 was 2358, checked in by josefine.ghattas, 10 years ago

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