1C Copyright 1981-2016 ECMWF. 2C 3C This software is licensed under the terms of the Apache Licence 4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5C 6C In applying this licence, ECMWF does not waive the privileges and immunities 7C granted to it by virtue of its status as an intergovernmental organisation 8C nor does it submit to any jurisdiction. 9C 10 11 INTEGER FUNCTION NEWISL(OLDGEO,NEWGEO,OLDLSM,OLDFLD,MASTER,NEWFLD) 12C 13C----> 14C**** NEWISL 15C 16C Purpose 17C ------- 18C 19C Interpolate a field based on an old land-sea mask to a field 20C based on a different land-sea mask. 21C 22C 23C Interface 24C --------- 25C 26C IRET = NEWISL(OLDGEO,NEWGEO,OLDLSM,OLDFLD,MASTER,NEWFLD) 27C 28C Input 29C ----- 30C 31C OLDGEO - GRIB section 2 describing grid of old field. 32C NEWGEO - GRIB section 2 describing grid of new field. 33C OLDLSM - Array of land-sea mask values for old field. 34C OLDFLD - Array of values for old field. 35C MASTER - Array of land-sea mask values for new field. 36C 37C 38C Output 39C ------ 40C 41C NEWFLD - Array of values for new field. 42C 43C Function returns: 44C - 0 if all is well 45C - 1, otherwise. 46C 47C 48C Method 49C ------ 50C 51C Build up field offsets from input geometries. 52C For each point of the new field, find in the old field the four 53C nearest neighbours' positions, values and types. 54C Calculate new value from the neighbours. 55C 56C 57C Externals 58C --------- 59C 60C IGGLAT - Compute gaussian lines of latitude. 61C JNORSGG - Find gaussian latitudes to north and south of latitude. 62C ISLPROC - Calculate value of new field point. 63C INTLOG - Log messages 64C 65C 66C Author 67C ------ 68C 69C J.D.Chambers ECMWF August 2000 70C 71C----< 72C 73 IMPLICIT NONE 74C 75C Function arguments 76C 77 INTEGER OLDGEO(*),NEWGEO(*) 78 REAL OLDLSM(*),OLDFLD(*),MASTER(*),NEWFLD(*) 79C 80#include "intisl.h" 81#include "parim.h" 82C 83C Parameters 84C 85 INTEGER JPMAXLT, JPGAUSS, JP1000 86 PARAMETER (JPMAXLT=721) 87 PARAMETER (JPGAUSS=4) 88 PARAMETER (JP1000=1000) 89C 90C Local variables 91C 92 INTEGER TOTAL, NEXT, LOOP, NEWOFF(JPMAXLT), OLDOFF(JPMAXLT) 93 INTEGER LATIT, LONG, NEWTYPE, IRET, NPTS 94 INTEGER LAT(2), LON(4) 95 INTEGER PT(4), TYPE(4) 96 REAL OLAT(2), OLON(4) 97 REAL RLATOLD(JPMAXLT), RLATNEW(JPMAXLT), RLAT, RLON 98C 99C Externals 100C 101 INTEGER IGGLAT, JNORSGG 102 REAL ISLPROC 103 EXTERNAL IGGLAT, JNORSGG, ISLPROC 104C 105C -----------------------------------------------------------------| 106C* Section 1. Build working values using input geometries. 107C -----------------------------------------------------------------| 108C 109 100 CONTINUE 110C 111 NEWISL = 0 112C 113C -----------------------------------------------------------------| 114C* Section 2. Calculate number of points in new field and offset 115C to start of each latitude in the new grid 116C -----------------------------------------------------------------| 117C 118 200 CONTINUE 119C 120 IF( NEWGEO(1).EQ.JPGAUSS ) THEN 121C 122C New field is gaussian 123C 124 CALL INTLOG(JP_DEBUG,'NEWISL: New field is gaussian',JPQUIET) 125C 126 IF( NEWGEO(17).EQ.0 ) THEN 127 CALL INTLOG(JP_DEBUG,'NEWISL: New field is regular',JPQUIET) 128 TOTAL = NEWGEO(2)*NEWGEO(3) 129 NEWOFF(1) = 0 130 DO LOOP = 2, NEWGEO(NJ) 131 NEWOFF(LOOP) = NEWOFF(LOOP-1) + NEWGEO(2) 132 ENDDO 133 ELSE 134 CALL INTLOG(JP_DEBUG,'NEWISL: New field is reduced',JPQUIET) 135 NEWOFF(1) = 0 136 TOTAL = NEWGEO(NPOINTS) 137 DO LOOP = 2, NEWGEO(NJ) 138 NEWOFF(LOOP) = NEWOFF(LOOP-1) + NEWGEO(NPOINTS-2+LOOP) 139 TOTAL = TOTAL + NEWGEO(NPOINTS-1+LOOP) 140 ENDDO 141 ENDIF 142C 143C Get the gaussian latitudes for the new field 144C 145 IRET = IGGLAT(NEWGEO(NGAUSS)*2,RLATNEW,0,-1) 146 IF( IRET.NE.0 ) THEN 147 WRITE(*,*) 'NEWISL: Problem call igglat for new grid' 148 NEWISL = 1 149 RETURN 150 ENDIF 151C 152 ELSE 153C 154C New field is lat/long 155C 156 CALL INTLOG(JP_DEBUG,'NEWISL: New field is lat/long',JPQUIET) 157 TOTAL = NEWGEO(2)*NEWGEO(3) 158 DO LOOP = 1, NEWGEO(3) 159 NEWOFF(LOOP) = NEWGEO(2)*(LOOP-1) 160 RLATNEW(LOOP) = 90.0 - (REAL((LOOP-1)*NEWGEO(10))/JP1000) 161 ENDDO 162C 163 ENDIF 164C 165 CALL INTLOG(JP_DEBUG,'NEWISL: No. of pts in new field = ',TOTAL) 166C 167C -----------------------------------------------------------------| 168C* Section 3. Get the gaussian latitudes for the old field and 169C setup the offsets to the start of each latitude. 170C -----------------------------------------------------------------| 171C 172 300 CONTINUE 173C 174 OLDOFF(1) = 0 175 DO LOOP = 2, OLDGEO(NJ) 176 OLDOFF(LOOP) = OLDOFF(LOOP-1) + OLDGEO(NPOINTS-2+LOOP) 177 ENDDO 178C 179 IRET = IGGLAT(OLDGEO(NGAUSS)*2,RLATOLD,0,-1) 180 IF( IRET.NE.0 ) THEN 181 WRITE(*,*) 'NEWISL: Problem call igglat for old grid' 182 NEWISL = 1 183 RETURN 184 ENDIF 185C 186C -----------------------------------------------------------------| 187C* Section 4. Work through the points in the new field. 188C -----------------------------------------------------------------| 189C 190 400 CONTINUE 191C 192 DO NEXT = 1, TOTAL 193C 194C Calculate lat/long 195C 196 DO LOOP = 1, NEWGEO(NJ) 197 IF( NEWOFF(LOOP).GE.NEXT ) THEN 198 LATIT = LOOP - 1 199 GOTO 410 200 ENDIF 201 ENDDO 202 LATIT = NEWGEO(NJ) 203C 204 410 CONTINUE 205C 206 LONG = NEXT - NEWOFF(LATIT) 207C 208 RLAT = RLATNEW(LATIT) 209 IF( NEWGEO(1).EQ.JPGAUSS ) THEN 210 IF( NEWGEO(17).EQ.0 ) THEN 211 RLON = REAL((LONG-1)*NEWGEO(9))/JP1000 212 ELSE 213 RLON = (REAL(LONG-1)*360.0)/REAL(NEWGEO(NPOINTS+LATIT-1)) 214 ENDIF 215 ELSE 216 RLON = REAL((LONG-1)*NEWGEO(9))/JP1000 217 ENDIF 218C 219C Find type of point (land or sea) 220C 221 IF( MASTER(NEXT).GT.MASTERTHRESHOLD ) THEN 222 NEWTYPE = LAND 223 ELSE 224 NEWTYPE = SEA 225 ENDIF 226C 227C Find four neighbours in the old field with their types 228C (Find NW neighbour and deduce the others). 229C 230 LAT(NORTH) = JNORSGG(RLAT,RLATOLD,OLDGEO(NGAUSS),1) 231 LAT(SOUTH) = JNORSGG(RLAT,RLATOLD,OLDGEO(NGAUSS),0) 232C 233 OLAT(NORTH) = RLATOLD(LAT(NORTH)) 234 OLAT(SOUTH) = RLATOLD(LAT(SOUTH)) 235C 236 NPTS = OLDGEO(NPOINTS-1+LAT(NORTH)) 237 LON(NWEST) = 1 + INT(RLON/(360.0/REAL(NPTS))) 238 LON(NEAST) = LON(NWEST) + 1 239 IF( LON(NEAST).GT.NPTS ) LON(NEAST) = 1 240C 241 OLON(NWEST) = (REAL(LON(NWEST)-1)*360.0)/REAL(NPTS) 242 IF( LON(NEAST).EQ.1 ) THEN 243 OLON(NEAST) = 360.0 244 ELSE 245 OLON(NEAST) = (REAL(LON(NEAST)-1)*360.0)/REAL(NPTS) 246 ENDIF 247C 248 NPTS = OLDGEO(NPOINTS-1+LAT(SOUTH)) 249 LON(SWEST) = 1 + INT(RLON/(360.0/REAL(NPTS))) 250 LON(SEAST) = LON(SWEST) + 1 251 IF( LON(SEAST).GT.NPTS ) LON(SEAST) = 1 252C 253 OLON(SWEST) = (REAL(LON(SWEST)-1)*360.0)/REAL(NPTS) 254 IF( LON(SEAST).EQ.1 ) THEN 255 OLON(SEAST) = 360.0 256 ELSE 257 OLON(SEAST) = (REAL(LON(SEAST)-1)*360.0)/REAL(NPTS) 258 ENDIF 259C 260 PT(NWEST) = OLDOFF(LAT(NORTH)) + LON(NWEST) 261 IF( OLDLSM(PT(NWEST)).GT.OLDLSMTHRESHOLD ) THEN 262 TYPE(NWEST) = LAND 263 ELSE 264 TYPE(NWEST) = SEA 265 ENDIF 266C 267 PT(NEAST) = OLDOFF(LAT(NORTH)) + LON(NEAST) 268 IF( OLDLSM(PT(NEAST)).GT.OLDLSMTHRESHOLD ) THEN 269 TYPE(NEAST) = LAND 270 ELSE 271 TYPE(NEAST) = SEA 272 ENDIF 273C 274 PT(SWEST) = OLDOFF(LAT(SOUTH)) + LON(SWEST) 275 IF( OLDLSM(PT(SWEST)).GT.OLDLSMTHRESHOLD ) THEN 276 TYPE(SWEST) = LAND 277 ELSE 278 TYPE(SWEST) = SEA 279 ENDIF 280C 281 PT(SEAST) = OLDOFF(LAT(SOUTH)) + LON(SEAST) 282 IF( OLDLSM(PT(SEAST)).GT.OLDLSMTHRESHOLD ) THEN 283 TYPE(SEAST) = LAND 284 ELSE 285 TYPE(SEAST) = SEA 286 ENDIF 287C 288C Interpolate the new value 289C 290 NEWFLD(NEXT) = 291 X ISLPROC(RLAT,RLON,NEWTYPE,OLAT,OLON,TYPE,PT,OLDFLD) 292C 293 ENDDO 294C 295C -----------------------------------------------------------------| 296C* Section 9. Return 297C -----------------------------------------------------------------| 298C 299 900 CONTINUE 300C 301C 302 RETURN 303 END 304 305