[7541] | 1 | MODULE interregxy |
---|
| 2 | |
---|
| 3 | ! Modules used : |
---|
| 4 | |
---|
| 5 | USE constantes |
---|
| 6 | USE mod_orchidee_para |
---|
| 7 | USE grid |
---|
| 8 | USE module_llxy |
---|
| 9 | USE polygones |
---|
| 10 | |
---|
| 11 | IMPLICIT NONE |
---|
| 12 | |
---|
| 13 | PRIVATE |
---|
| 14 | PUBLIC interregxy_aggr2d, interregxy_aggrve |
---|
| 15 | |
---|
| 16 | CONTAINS |
---|
| 17 | |
---|
| 18 | SUBROUTINE interregxy_aggr2d(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
| 19 | & iml, jml, lon_rel, lat_rel, mask, callsign, & |
---|
| 20 | & incmax, indinc, areaoverlap, ok) |
---|
| 21 | ! |
---|
| 22 | ! INPUT |
---|
| 23 | ! |
---|
| 24 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
| 25 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
| 26 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
| 27 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
| 28 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
| 29 | INTEGER(i_std), INTENT(in) :: iml, jml ! Size of the source grid |
---|
| 30 | REAL(r_std), INTENT(in) :: lon_rel(iml, jml) ! Longitudes for the source grid |
---|
| 31 | REAL(r_std), INTENT(in) :: lat_rel(iml, jml) ! Latitudes for the souce grid |
---|
| 32 | INTEGER(i_std), INTENT(in) :: mask(iml, jml) ! Mask which retains only the significative points |
---|
| 33 | ! of the source grid. |
---|
| 34 | CHARACTER(LEN=*), INTENT(in) :: callsign ! Allows to specify which variable is beeing treated |
---|
| 35 | INTEGER(i_std), INTENT(in) :: incmax ! Maximum point of the source grid we can store. |
---|
| 36 | ! |
---|
| 37 | ! Output |
---|
| 38 | ! |
---|
| 39 | INTEGER(i_std), INTENT(out) :: indinc(nbpt,incmax,2) |
---|
| 40 | REAL(r_std), INTENT(out) :: areaoverlap(nbpt,incmax) |
---|
| 41 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
| 42 | ! |
---|
| 43 | ! Local Variables |
---|
| 44 | ! |
---|
| 45 | ! The number of sub-division we wish to have on each side of the rectangle. |
---|
| 46 | ! sbd=1 means that on each side we take only the middle point. |
---|
| 47 | ! |
---|
| 48 | INTEGER(i_std), PARAMETER :: sbd=5 |
---|
| 49 | INTEGER(i_std), PARAMETER :: idom=1 |
---|
| 50 | ! |
---|
| 51 | INTEGER(i_std) :: i, is, js, in |
---|
| 52 | INTEGER(i_std) :: ilow, iup, jlow, jup |
---|
| 53 | REAL(r_std) :: dle, dlw, dls, dln |
---|
| 54 | ! |
---|
| 55 | REAL(r_std), DIMENSION((sbd+1)*4+1) :: lons, lats, ri, rj |
---|
| 56 | ! |
---|
| 57 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: nbov |
---|
| 58 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sourcei, sourcej |
---|
| 59 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: overlap |
---|
| 60 | LOGICAL :: checkprint=.FALSE. |
---|
| 61 | ! |
---|
| 62 | ! |
---|
| 63 | REAL(r_std), DIMENSION(2) :: BBlon, BBlat |
---|
| 64 | ! |
---|
| 65 | WRITE(*,*) "interregxy_aggr2d working on ", callsign |
---|
| 66 | ! |
---|
| 67 | ALLOCATE(sourcei(iim_g, jjm_g, incmax)) |
---|
| 68 | ALLOCATE(sourcej(iim_g, jjm_g, incmax)) |
---|
| 69 | ALLOCATE(nbov(iim_g, jjm_g)) |
---|
| 70 | ALLOCATE(overlap(iim_g, jjm_g, incmax)) |
---|
| 71 | nbov(:,:) = 0 |
---|
| 72 | ! |
---|
| 73 | BBlon(1) = MINVAL(lalo(:,2)) |
---|
| 74 | BBlon(2) = MAXVAL(lalo(:,2)) |
---|
| 75 | BBlat(1) = MINVAL(lalo(:,1)) |
---|
| 76 | BBlat(2) = MAXVAL(lalo(:,1)) |
---|
| 77 | ! |
---|
| 78 | ! Go through the entire source grid and try to place each grid box in a target grid box. |
---|
| 79 | ! |
---|
| 80 | DO is=1,iml |
---|
| 81 | DO js=1,jml |
---|
| 82 | ! |
---|
| 83 | ! 1.0 Get the center of the box and their corresponding grid point |
---|
| 84 | ! |
---|
| 85 | lons(1) = lon_rel(is,js) |
---|
| 86 | lats(1) = lat_rel(is,js) |
---|
| 87 | ! |
---|
| 88 | ! Assume all grids go from West to Est |
---|
| 89 | dle=(lon_rel(MIN(is+1,iml),js)-lon_rel(is,js))/2.0 |
---|
| 90 | dlw=(lon_rel(is,js)-lon_rel(MAX(is-1,1),js))/2.0 |
---|
| 91 | ! |
---|
| 92 | IF ( lat_rel(1,1) < lat_rel(iml,jml)) THEN |
---|
| 93 | ! Grid in SN direction |
---|
| 94 | dls=(lat_rel(is,js)-lat_rel(is,MAX(js-1,1)))/2.0 |
---|
| 95 | dln=(lat_rel(is,MIN(js+1,jml))-lat_rel(is,js))/2.0 |
---|
| 96 | ELSE |
---|
| 97 | dls=(lat_rel(is,js)-lat_rel(is,MIN(js+1,jml)))/2.0 |
---|
| 98 | dln=(lat_rel(is,MAX(js-1,1))-lat_rel(is,js))/2.0 |
---|
| 99 | ENDIF |
---|
| 100 | ! Treat the border cases |
---|
| 101 | IF ( dlw < dle/2. ) dlw = dle |
---|
| 102 | IF ( dle < dlw/2. ) dle = dlw |
---|
| 103 | ! |
---|
| 104 | IF ( dls < dln/2. ) dls = dln |
---|
| 105 | IF ( dln < dls/2. ) dln = dls |
---|
| 106 | ! |
---|
| 107 | DO i=0,sbd |
---|
| 108 | ! Put the points in the order : SW, SE, NE, NW |
---|
| 109 | lons(2+i) = lons(1) - dlw + i*((dlw+dle)/FLOAT(sbd+1)) |
---|
| 110 | lats(2+i) = lats(1) - dls |
---|
| 111 | ! |
---|
| 112 | lons(2+(sbd+1)+i) = lons(1) + dle |
---|
| 113 | lats(2+(sbd+1)+i) = lats(1) - dls + i*((dls+dln)/FLOAT(sbd+1)) |
---|
| 114 | ! |
---|
| 115 | lons(2+2*(sbd+1)+i) = lons(1) + dle - i*((dlw+dle)/FLOAT(sbd+1)) |
---|
| 116 | lats(2+2*(sbd+1)+i) = lats(1) + dln |
---|
| 117 | ! |
---|
| 118 | lons(2+3*(sbd+1)+i) = lons(1) - dlw |
---|
| 119 | lats(2+3*(sbd+1)+i) = lats(1) + dln - i*((dls+dln)/FLOAT(sbd+1)) |
---|
| 120 | ENDDO |
---|
| 121 | ! |
---|
| 122 | ! Test that the polygone overlaps with the Bounding box. We might have projections |
---|
| 123 | ! which are not valid outside of the RCM domain. |
---|
| 124 | ! |
---|
| 125 | IF ( ((MINVAL(lons) > BBlon(1) .AND. MINVAL(lons) < BBlon(2)) .OR. & |
---|
| 126 | & (MAXVAL(lons) > BBlon(1) .AND. MAXVAL(lons) < BBlon(2))) .AND. & |
---|
| 127 | & ((MINVAL(lats) > BBlat(1) .AND. MINVAL(lats) < BBlat(2)) .OR. & |
---|
| 128 | & (MAXVAL(lats) > BBlat(1) .AND. MAXVAL(lats) < BBlat(2)))) THEN |
---|
| 129 | ! |
---|
| 130 | ! |
---|
| 131 | CALL Grid_Toij (lons, lats, ri, rj) |
---|
| 132 | ! |
---|
| 133 | ! |
---|
| 134 | ! Find the points of the WRF grid boxes we need to scan to place this box (is, js) |
---|
| 135 | ! of the source grid. |
---|
| 136 | ! |
---|
| 137 | ilow = MAX(FLOOR(MINVAL(ri)),1) |
---|
| 138 | iup = MIN(CEILING(MAXVAL(ri)), iim_g) |
---|
| 139 | jlow = MAX(FLOOR(MINVAL(rj)),1) |
---|
| 140 | jup = MIN(CEILING(MAXVAL(rj)), jjm_g) |
---|
| 141 | ! |
---|
| 142 | ! |
---|
| 143 | IF ( ilow < iup .OR. jlow < jup ) THEN |
---|
| 144 | ! |
---|
| 145 | ! Find the grid boxes in WRF and compute the overlap area. |
---|
| 146 | ! |
---|
| 147 | CALL interregxy_findpoints(iim_g, jjm_g, lon_rel(is,js), lon_rel(is,js), is, js, ri, rj, incmax, & |
---|
| 148 | & checkprint, nbov, sourcei, sourcej, overlap) |
---|
| 149 | |
---|
| 150 | ENDIF |
---|
| 151 | ! |
---|
| 152 | ! |
---|
| 153 | ENDIF |
---|
| 154 | ENDDO |
---|
| 155 | ENDDO |
---|
| 156 | ! |
---|
| 157 | ! Gather the land points from the global grid |
---|
| 158 | ! |
---|
| 159 | areaoverlap(:,:) = zero |
---|
| 160 | DO in=1,nbpt |
---|
| 161 | ! |
---|
| 162 | is = ilandindex(in) |
---|
| 163 | js = jlandindex(in) |
---|
| 164 | ! |
---|
| 165 | indinc(in,1:nbov(is,js),1) = sourcei(is,js,1:nbov(is,js)) |
---|
| 166 | indinc(in,1:nbov(is,js),2) = sourcej(is,js,1:nbov(is,js)) |
---|
| 167 | areaoverlap(in,1:nbov(is,js)) = overlap(is,js,1:nbov(is,js)) |
---|
| 168 | ENDDO |
---|
| 169 | ! |
---|
| 170 | ! The model will stop if incmax is not sufficient. So if we are here all is OK. |
---|
| 171 | ! |
---|
| 172 | ok=.TRUE. |
---|
| 173 | ! |
---|
| 174 | END SUBROUTINE interregxy_aggr2d |
---|
| 175 | ! |
---|
| 176 | ! =================================================================================== |
---|
| 177 | ! |
---|
| 178 | SUBROUTINE interregxy_aggrve(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
| 179 | & iml, lon_rel, lat_rel, resol_lon, resol_lat, callsign, & |
---|
| 180 | & incmax, indinc, areaoverlap, ok) |
---|
| 181 | ! |
---|
| 182 | ! INPUT |
---|
| 183 | ! |
---|
| 184 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
| 185 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
| 186 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
| 187 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
| 188 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
| 189 | INTEGER(i_std), INTENT(in) :: iml ! Size of the source grid |
---|
| 190 | REAL(r_std), INTENT(in) :: lon_rel(iml) ! Longitudes for the source grid |
---|
| 191 | REAL(r_std), INTENT(in) :: lat_rel(iml) ! Latitudes for the source grid |
---|
| 192 | REAL(r_std), INTENT(in) :: resol_lon, resol_lat ! Resolution in meters of the source grid |
---|
| 193 | CHARACTER(LEN=*), INTENT(in) :: callsign ! Allows to specify which variable is beeing treated |
---|
| 194 | INTEGER(i_std), INTENT(in) :: incmax ! Maximum point of the source grid we can store. |
---|
| 195 | ! |
---|
| 196 | ! Output |
---|
| 197 | ! |
---|
| 198 | INTEGER(i_std), INTENT(out) :: indinc(nbpt,incmax) |
---|
| 199 | REAL(r_std), INTENT(out) :: areaoverlap(nbpt,incmax) |
---|
| 200 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
| 201 | ! |
---|
| 202 | ! |
---|
| 203 | ! Local Variables |
---|
| 204 | ! |
---|
| 205 | ! The number of sub-division we wish to have on each side of the rectangle. |
---|
| 206 | ! sbd=1 means that on each side we take only the middle point. |
---|
| 207 | ! |
---|
| 208 | INTEGER(i_std), PARAMETER :: sbd=5 |
---|
| 209 | INTEGER(i_std), PARAMETER :: idom=1 |
---|
| 210 | ! |
---|
| 211 | INTEGER(i_std) :: i, is, js, in |
---|
| 212 | INTEGER(i_std) :: ilow, iup, jlow, jup |
---|
| 213 | REAL(r_std) :: dle, dlw, dls, dln, coslat |
---|
| 214 | ! |
---|
| 215 | REAL(r_std), DIMENSION((sbd+1)*4+1) :: lons, lats, ri, rj |
---|
| 216 | ! |
---|
| 217 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: nbov |
---|
| 218 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sourcei, sourcej |
---|
| 219 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: overlap |
---|
| 220 | ! |
---|
| 221 | ! |
---|
| 222 | REAL(r_std), DIMENSION(2) :: BBlon, BBlat |
---|
| 223 | REAL(r_std) :: half_resol_lon, half_resol_lat, trtmp |
---|
| 224 | LOGICAL :: checkprint=.FALSE. |
---|
| 225 | ! |
---|
| 226 | WRITE(*,*) "interregxy_aggrve working on ", callsign |
---|
| 227 | ! |
---|
| 228 | ALLOCATE(sourcei(iim_g, jjm_g, incmax)) |
---|
| 229 | ALLOCATE(sourcej(iim_g, jjm_g, incmax)) |
---|
| 230 | ALLOCATE(nbov(iim_g, jjm_g)) |
---|
| 231 | ALLOCATE(overlap(iim_g, jjm_g, incmax)) |
---|
| 232 | nbov(:,:) = 0 |
---|
| 233 | ! |
---|
| 234 | BBlon(1) = MINVAL(lalo(:,2)) |
---|
| 235 | BBlon(2) = MAXVAL(lalo(:,2)) |
---|
| 236 | BBlat(1) = MINVAL(lalo(:,1)) |
---|
| 237 | BBlat(2) = MAXVAL(lalo(:,1)) |
---|
| 238 | ! |
---|
| 239 | ! Go through the entire source grid and try to place each grid box in a target grid box. |
---|
| 240 | ! |
---|
| 241 | DO is=1,iml |
---|
| 242 | ! |
---|
| 243 | ! 1.0 Get the center of the box and their corresponding grid point |
---|
| 244 | ! |
---|
| 245 | lons(1) = lon_rel(is) |
---|
| 246 | lats(1) = lat_rel(is) |
---|
| 247 | ! Some pre-calculations |
---|
| 248 | trtmp = 180./(2.*pi) |
---|
| 249 | half_resol_lon = resol_lon/2.0 |
---|
| 250 | half_resol_lat = resol_lat/2.0 |
---|
| 251 | ! |
---|
| 252 | dln = half_resol_lat*trtmp/R_Earth |
---|
| 253 | dls = half_resol_lat*trtmp/R_Earth |
---|
| 254 | ! |
---|
| 255 | DO i=0,sbd |
---|
| 256 | ! Put the points in the order : SW, SE, NE, NW |
---|
| 257 | lats(2+i) = lats(1) - dls |
---|
| 258 | coslat = COS( lats(2+i) * pi/180. ) |
---|
| 259 | dlw = half_resol_lon*trtmp/R_Earth/coslat |
---|
| 260 | lons(2+i) = lons(1) - dlw + i*((2.*dlw)/FLOAT(sbd+1)) |
---|
| 261 | ! |
---|
| 262 | lats(2+(sbd+1)+i) = lats(1) - dls + i*((dls+dln)/FLOAT(sbd+1)) |
---|
| 263 | coslat = COS( lats(2+(sbd+1)+i) * pi/180. ) |
---|
| 264 | dle = half_resol_lon*trtmp/R_Earth/coslat |
---|
| 265 | lons(2+(sbd+1)+i) = lons(1) + dle |
---|
| 266 | ! |
---|
| 267 | lats(2+2*(sbd+1)+i) = lats(1) + dln |
---|
| 268 | coslat = COS( lats(2+2*(sbd+1)+i) * pi/180. ) |
---|
| 269 | dle = half_resol_lon*trtmp/R_Earth/coslat |
---|
| 270 | lons(2+2*(sbd+1)+i) = lons(1) + dle - i*((2*dle)/FLOAT(sbd+1)) |
---|
| 271 | ! |
---|
| 272 | lats(2+3*(sbd+1)+i) = lats(1) + dln - i*((dls+dln)/FLOAT(sbd+1)) |
---|
| 273 | coslat = COS( lats(2+3*(sbd+1)+i) * pi/180. ) |
---|
| 274 | dlw = half_resol_lon*trtmp/R_Earth/coslat |
---|
| 275 | lons(2+3*(sbd+1)+i) = lons(1) - dlw |
---|
| 276 | ENDDO |
---|
| 277 | ! |
---|
| 278 | ! Test that the polygone overlaps with the Bounding box. We might have projections |
---|
| 279 | ! which are not valid outside of the RCM domain. |
---|
| 280 | ! |
---|
| 281 | IF ( ((MINVAL(lons) > BBlon(1) .AND. MINVAL(lons) < BBlon(2)) .OR. & |
---|
| 282 | & (MAXVAL(lons) > BBlon(1) .AND. MAXVAL(lons) < BBlon(2))) .AND. & |
---|
| 283 | & ((MINVAL(lats) > BBlat(1) .AND. MINVAL(lats) < BBlat(2)) .OR. & |
---|
| 284 | & (MAXVAL(lats) > BBlat(1) .AND. MAXVAL(lats) < BBlat(2)))) THEN |
---|
| 285 | ! |
---|
| 286 | CALL Grid_Toij (lons, lats, ri, rj) |
---|
| 287 | ! |
---|
| 288 | ! |
---|
| 289 | ! Find the points of the WRF grid boxes we need to scan to place this box (is, js) |
---|
| 290 | ! of the source grid. |
---|
| 291 | ! |
---|
| 292 | ilow = MAX(FLOOR(MINVAL(ri)),1) |
---|
| 293 | iup = MIN(CEILING(MAXVAL(ri)), iim_g) |
---|
| 294 | jlow = MAX(FLOOR(MINVAL(rj)),1) |
---|
| 295 | jup = MIN(CEILING(MAXVAL(rj)), jjm_g) |
---|
| 296 | ! |
---|
| 297 | ! |
---|
| 298 | IF ( ilow < iup .OR. jlow < jup ) THEN |
---|
| 299 | ! |
---|
| 300 | ! Find the grid boxes in WRF and compute the overlap area. |
---|
| 301 | ! |
---|
| 302 | CALL interregxy_findpoints(iim_g, jjm_g, lon_rel(is), lon_rel(is), is, 1, ri, rj, incmax, & |
---|
| 303 | & checkprint, nbov, sourcei, sourcej, overlap) |
---|
| 304 | |
---|
| 305 | ENDIF |
---|
| 306 | ! |
---|
| 307 | ! |
---|
| 308 | ENDIF |
---|
| 309 | ENDDO |
---|
| 310 | ! |
---|
| 311 | ! Gather the land points from the global grid |
---|
| 312 | ! |
---|
| 313 | areaoverlap(:,:) = zero |
---|
| 314 | DO in=1,nbpt |
---|
| 315 | ! |
---|
| 316 | is = ilandindex(in) |
---|
| 317 | js = jlandindex(in) |
---|
| 318 | ! |
---|
| 319 | indinc(in,1:nbov(is,js)) = sourcei(is,js,1:nbov(is,js)) |
---|
| 320 | areaoverlap(in,1:nbov(is,js)) = overlap(is,js,1:nbov(is,js)) |
---|
| 321 | ENDDO |
---|
| 322 | ! |
---|
| 323 | ! The model will stop if incmax is not sufficient. So if we are here all is OK. |
---|
| 324 | ! |
---|
| 325 | ok=.TRUE. |
---|
| 326 | ! |
---|
| 327 | END SUBROUTINE interregxy_aggrve |
---|
| 328 | ! |
---|
| 329 | ! =================================================================================== |
---|
| 330 | ! |
---|
| 331 | SUBROUTINE interregxy_findpoints(nb_we, nb_sn, lon_cen, lat_cen, is, js, ri, rj, & |
---|
| 332 | & nbpt_max, checkprint, nbpt, sourcei, sourcej, overlap_area) |
---|
| 333 | ! |
---|
| 334 | ! |
---|
| 335 | ! This Subroutine finds the overlap of grid box ri,rj with the WRF grid. (ri,rj) are the WRF indexes of the (is,js) point |
---|
| 336 | ! of the source grid. |
---|
| 337 | ! Each time we find an overlap of (is,js) with a WRF grid box (ix,jx) we add that in the (sourcei,sourcej) fields. |
---|
| 338 | ! |
---|
| 339 | ! |
---|
| 340 | ! Arguments |
---|
| 341 | ! |
---|
| 342 | INTEGER(i_std), INTENT(in) :: nb_we, nb_sn |
---|
| 343 | REAL(r_std), INTENT(in) :: lon_cen, lat_cen |
---|
| 344 | INTEGER(i_std), INTENT(in) :: is, js |
---|
| 345 | REAL(r_std), DIMENSION(:), INTENT(in) :: ri, rj |
---|
| 346 | INTEGER(i_std), INTENT(in) :: nbpt_max |
---|
| 347 | LOGICAL, INTENT(in) :: checkprint |
---|
| 348 | INTEGER(i_std), DIMENSION(:,:), INTENT(out) :: nbpt |
---|
| 349 | INTEGER(i_std), DIMENSION(:,:,:), INTENT(out) :: sourcei, sourcej |
---|
| 350 | REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: overlap_area |
---|
| 351 | ! |
---|
| 352 | ! Local |
---|
| 353 | ! |
---|
| 354 | INTEGER(i_std) :: ilow, iup, jlow, jup, nbr |
---|
| 355 | INTEGER(i_std) :: ix, jx, i |
---|
| 356 | |
---|
| 357 | INTEGER(i_std) :: nbint, nbcross, nbsorted, nbfinal, nbtmpa, nbtmpb |
---|
| 358 | ! |
---|
| 359 | INTEGER(i_std) :: nbdots=5, dimpoly = 200 |
---|
| 360 | REAL(r_std), DIMENSION(4,2) :: poly_wrf |
---|
| 361 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: poly_in, poly_tmpa, poly_tmpb, poly_int, poly_cross |
---|
| 362 | ! |
---|
| 363 | REAL(r_std) :: ax, ay, overlaparea |
---|
| 364 | ! |
---|
| 365 | INTEGER(i_std), PARAMETER :: icheck=50, jcheck=46 |
---|
| 366 | ! |
---|
| 367 | ! |
---|
| 368 | ! |
---|
| 369 | nbr = SIZE(ri) |
---|
| 370 | ! |
---|
| 371 | IF ( dimpoly < 4*(nbr-1) ) THEN |
---|
| 372 | dimpoly = 4*nbr |
---|
| 373 | ENDIF |
---|
| 374 | ! |
---|
| 375 | ALLOCATE(poly_in(nbr-1,2)) |
---|
| 376 | ALLOCATE(poly_tmpa(dimpoly,2)) |
---|
| 377 | ALLOCATE(poly_tmpb(dimpoly,2)) |
---|
| 378 | ALLOCATE(poly_int(dimpoly,2)) |
---|
| 379 | ALLOCATE(poly_cross(dimpoly,2)) |
---|
| 380 | ! |
---|
| 381 | ! Scan all Determine the range of indicis of the WRF grid we need to scan |
---|
| 382 | ! within the box of the source grid. |
---|
| 383 | ! |
---|
| 384 | ilow = MAX(FLOOR(MINVAL(ri)),1) |
---|
| 385 | iup = MIN(CEILING(MAXVAL(ri)), nb_we) |
---|
| 386 | jlow = MAX(FLOOR(MINVAL(rj)),1) |
---|
| 387 | jup = MIN(CEILING(MAXVAL(rj)), nb_sn) |
---|
| 388 | ! |
---|
| 389 | ! |
---|
| 390 | poly_wrf(:,1) = (/ -0.5, 0.5, 0.5, -0.5 /) |
---|
| 391 | poly_wrf(:,2) = (/ -0.5, -0.5, 0.5, 0.5 /) |
---|
| 392 | ! |
---|
| 393 | DO ix = ilow,iup |
---|
| 394 | DO jx = jlow,jup |
---|
| 395 | ! |
---|
| 396 | ! Bring the points of the source grid into the +-0.5 range of the WRF grid ... |
---|
| 397 | ! and remove the center point. |
---|
| 398 | ! |
---|
| 399 | poly_in(1:nbr-1,1) = ri(2:nbr) - ix |
---|
| 400 | poly_in(1:nbr-1,2) = rj(2:nbr) - jx |
---|
| 401 | ! |
---|
| 402 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 403 | WRITE(*,*) "X = ", poly_in(1:nbr-1,1) |
---|
| 404 | WRITE(*,*) "Y = ", poly_in(1:nbr-1,2) |
---|
| 405 | ENDIF |
---|
| 406 | ! |
---|
| 407 | CALL polygones_extend(nbr-1, poly_in, nbdots, nbtmpa, poly_tmpa) |
---|
| 408 | CALL polygones_extend(4, poly_wrf, nbdots, nbtmpb, poly_tmpb) |
---|
| 409 | ! |
---|
| 410 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 411 | WRITE(*,*) "X_extend = ", poly_tmpa(1:nbtmpa,1) |
---|
| 412 | WRITE(*,*) "Y_extend = ", poly_tmpa(1:nbtmpa,2) |
---|
| 413 | ENDIF |
---|
| 414 | ! |
---|
| 415 | CALL polygones_intersection(nbtmpa, poly_tmpa, nbtmpb, poly_tmpb, nbint, poly_int) |
---|
| 416 | CALL polygones_crossing(nbr-1, poly_in, 4, poly_wrf, nbcross, poly_cross) |
---|
| 417 | ! |
---|
| 418 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 419 | WRITE(*,*) "nbint = nbcross = ", nbint, nbcross |
---|
| 420 | ENDIF |
---|
| 421 | ! |
---|
| 422 | IF ( nbint+nbcross > dimpoly ) THEN |
---|
| 423 | WRITE(*,*) "The intersection of polygones is larger than the memory foressen : ", nbint+nbcross, dimpoly |
---|
| 424 | STOP "interregxy_compoverlap" |
---|
| 425 | ENDIF |
---|
| 426 | ! |
---|
| 427 | ! |
---|
| 428 | IF ( nbint+nbcross > 2 ) THEN |
---|
| 429 | |
---|
| 430 | poly_tmpa(1:nbint,:) = poly_int(1:nbint,:) |
---|
| 431 | poly_tmpa(nbint+1:nbint+nbcross,:) = poly_cross(1:nbcross,:) |
---|
| 432 | |
---|
| 433 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 434 | WRITE(*,*) "X_overlap = ", poly_tmpa(1:(nbint+nbcross),1) |
---|
| 435 | WRITE(*,*) "Y_overlap = ", poly_tmpa(1:(nbint+nbcross),2) |
---|
| 436 | ENDIF |
---|
| 437 | |
---|
| 438 | CALL polygones_cleanup(nbint+nbcross, poly_tmpa, nbfinal, poly_tmpb) |
---|
| 439 | IF ( nbint+nbcross < 3 ) THEN |
---|
| 440 | WRITE(*,*) "poly_tmpa X : ", poly_tmpa(1:(nbint+nbcross),1) |
---|
| 441 | WRITE(*,*) "poly_tmpa Y : ", poly_tmpa(1:(nbint+nbcross),2) |
---|
| 442 | WRITE(*,*) "Cleaning goes too far : ", nbint+nbcross, " to ", nbfinal |
---|
| 443 | WRITE(*,*) "poly_tmpb X : ", poly_tmpb(1:nbfinal,1) |
---|
| 444 | WRITE(*,*) "poly_tmpb Y : ", poly_tmpb(1:nbfinal,2) |
---|
| 445 | ENDIF |
---|
| 446 | |
---|
| 447 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 448 | WRITE(*,*) "X_final = ", poly_tmpb(1:nbfinal,1) |
---|
| 449 | WRITE(*,*) "Y_final = ", poly_tmpb(1:nbfinal,2) |
---|
| 450 | ENDIF |
---|
| 451 | IF ( nbfinal > 2 ) THEN |
---|
| 452 | CALL polygones_convexhull(nbfinal, poly_tmpb, nbsorted, poly_tmpa) |
---|
| 453 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 454 | WRITE(*,*) "X_sorted = ", poly_tmpa(1:nbsorted,1) |
---|
| 455 | WRITE(*,*) "Y_sorted = ", poly_tmpa(1:nbsorted,2) |
---|
| 456 | ENDIF |
---|
| 457 | CALL polygones_area(nbsorted, poly_tmpa, dxWRF(ix,jx), dyWRF(ix,jx), overlaparea) |
---|
| 458 | ELSE |
---|
| 459 | overlaparea = 0.0 |
---|
| 460 | ENDIF |
---|
| 461 | |
---|
| 462 | IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN |
---|
| 463 | WRITE(*,*) ix, jx, "overlaparea = ", overlaparea |
---|
| 464 | WRITE(*,*) "=============================================================================" |
---|
| 465 | ENDIF |
---|
| 466 | |
---|
| 467 | IF ( overlaparea > 0.0 ) THEN |
---|
| 468 | nbpt(ix,jx) = nbpt(ix,jx)+1 |
---|
| 469 | IF ( nbpt(ix,jx) > nbpt_max ) THEN |
---|
| 470 | WRITE(*,*) "ERROR in interregxy_findpoints. " |
---|
| 471 | WRITE(*,*) "Consider your algorithm to estimate nbpt_max. " |
---|
| 472 | WRITE(*,*) "You should change it as it does not give enough points." |
---|
| 473 | WRITE(*,*) "nbpt(ix,jx) = ", nbpt(ix,jx) |
---|
| 474 | WRITE(*,*) "nbpt_max = ", nbpt_max |
---|
| 475 | STOP "interregxy_findpoints" |
---|
| 476 | ENDIF |
---|
| 477 | sourcei(ix,jx,nbpt(ix,jx)) = is |
---|
| 478 | sourcej(ix,jx,nbpt(ix,jx)) = js |
---|
| 479 | overlap_area(ix,jx,nbpt(ix,jx)) = overlaparea |
---|
| 480 | ! |
---|
| 481 | ENDIF |
---|
| 482 | ENDIF |
---|
| 483 | ENDDO |
---|
| 484 | ENDDO |
---|
| 485 | ! |
---|
| 486 | DEALLOCATE(poly_in, poly_tmpa, poly_tmpb, poly_int, poly_cross) |
---|
| 487 | ! |
---|
| 488 | END SUBROUTINE interregxy_findpoints |
---|
| 489 | |
---|
| 490 | END MODULE interregxy |
---|