source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_global/interpol_help.f90 @ 7346

Last change on this file since 7346 was 5667, checked in by anne-sofie.lanso, 6 years ago

DEV: This revision allows for the usage of the new N input files used by ORCHIDEE-CN. In order to do so major revisions has been made to the code in slowproc in particual slowproc_ninput, slowproc_init and slowproc_main. Minor changes have also been made in nitrogen_dynamics. Moreover, a problem related to XIOS and the lack of a defined array for outputs with the dimension of the short product pool, nshort, was fixed. The problem arose because nshort was changed from 1 to 2 in r5656. In order to run this revision, you need to adjust your configuration ton the configeration for ORCHIDEE-CN-CAN found on the wiki: ​https://forge.ipsl.jussieu.fr/orchidee/browser/branches/ORCHIDEE-CN-CAN/ORCHIDEE_OL.

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