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