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