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 WV2DINT(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE, 12 X NEWWAVE,NORTH,WEST,KNSPEC,PMISS,RNS) 13C 14C----> 15C*****WV2DINT* 16C 17C PURPOSE 18C ------- 19C 20C To interpolate a wave field quasi-regular latitude 21C longitude grid to a regular latitude/longitude grid. 22C 23C 24C INTERFACE 25C --------- 26C 27C IRET = WV2DINT(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE, 28C NEWWAVE,NORTH,WEST,KNSPEC,PMISS,RNS) 29C 30C Input arguments 31C --------------- 32C 33C KNUM - No. of meridians from North to South pole (input field) 34C NUMPTS - Array giving number of points along each latitude 35C (empty latitudes have entry 0) 36C KE_W - First dimension of new array 37C = Number of points E-W in new grid 38C KN_S - Second dimension of new array 39C = Number of points N-S in new grid 40C RESON - Output grid resolution (degrees) 41C OLDWAVE - Original wave field 42C NORTH - Output grid northernmost latitude (degree) 43C WEST - Output grid westernmost longitude (degree) 44C KNSPEC - Number of 2D spectra values at each wave point 45C PMISS - Missing data value 46C RNS - Difference in degrees in NS disrection 47C 48C Output arguments 49C ---------------- 50 51C NEWWAVE - New wave field 52C 53C Function returns 0 if the interpolation was OK. 54C 55C 56C METHOD 57C ------ 58C 59C The nearest neighbouring grid point value is assigned to 60C the interpolated point. 61C 62C The output (sub) is the same as the input area. 63C 64C 65C EXTERNALS 66C --------- 67C 68C INTLOG - Log error message. 69C 70C 71C REFERENCE 72C --------- 73C 74C Based on: 75C SUBROUTINE INTERPOLATE 76C Peter Janssen ECMWF September 1995 77C and: 78C SUBROUTINE EXPOINT 79C Heinz Gunther ECMWF December 1989 80C 81C 82C Author. 83C ------- 84C 85C J.D.Chambers ECMWF November 1996 86C 87C J.D.Chambers ECMWF September 1998 88C Modified to handle subarea input/output. 89C The output (sub) is the same as the input area. 90C 91C 92C----< 93C 94 IMPLICIT NONE 95C 96#include "parim.h" 97C 98C Parameters 99 INTEGER JPROUTINE 100 PARAMETER ( JPROUTINE = 19410 ) 101 INTEGER JPLLMAX 102 PARAMETER ( JPLLMAX = 1801 ) 103C `--> allow upto 0.1 degree resolution for 104C Mediterranean, 0.1 for Global 105 INTEGER JPNMOUT 106 PARAMETER ( JPNMOUT = 1800 ) 107C `--> allow upto 0.1 degree resolution for 108C Mediterranean, 0.1 for Global 109C 110C Subroutine arguments 111C 112 INTEGER KNUM, NUMPTS, KE_W, KN_S, KNSPEC 113 DIMENSION NUMPTS(*) 114 REAL RESON, OLDWAVE, NEWWAVE, NORTH, WEST, PMISS 115 REAL RNS 116 DIMENSION OLDWAVE(*) 117 DIMENSION NEWWAVE(KE_W,KN_S) 118C 119C Local arguments 120C 121 INTEGER NLAT, K, NEWCOL, NEWROW 122 INTEGER I_N, I_S 123 REAL DELAT 124 REAL DIST_N, DIST_S 125 REAL LAT 126 DIMENSION LAT(JPNMOUT+1) 127 INTEGER INDEX 128 DIMENSION INDEX(JPLLMAX) 129C 130 REAL OLDLATS(JPLLMAX), OLDEAST 131 INTEGER NEXT, NUMNEW(JPLLMAX) 132 INTEGER IRET, W251IDX 133 EXTERNAL W251IDX 134 INTEGER XKNUM,XKE_W,XKN_S 135 DATA XKNUM/-1/,XKE_W/-1/,XKN_S/-1/ 136 REAL XRESON,XNORTH,XWEST 137 DATA XRESON/-999.0/,XNORTH/-999.0/,XWEST/-999.0/ 138 INTEGER IFIRST, NEWIDX(3600*1801) 139 DATA IFIRST/0/ 140 REAL DEPS 141 DATA DEPS/0.00005/ 142C 143 SAVE IFIRST,XKNUM,XKE_W,XKN_S,XRESON,XNORTH,XWEST,NEWIDX 144C 145C 146C --------------------------------------------------------------------- 147C* Section 1. Initalisation. 148C --------------------------------------------------------------------- 149C 150 100 CONTINUE 151C 152 WV2DINT = 0 153 CALL INTLOG(JP_DEBUG, 154 X 'WV2DINT: Wave interpolation requested.',JPQUIET) 155C 156C Check only new-style 2D wave spectra (parameter 251) 157C 158 IF( KNSPEC.GT.1 ) THEN 159 CALL INTLOG(JP_ERROR, 160 X 'WV2DINT: Only single-value 2D spectra field allowed',JPQUIET) 161 CALL INTLOG(JP_ERROR, 162 X 'WV2DINT: Value given = ',KNSPEC) 163 WV2DINT = JPROUTINE + 1 164 GOTO 900 165 ENDIF 166C 167C Check reduced latitude/longitude grid specification 168C 169 IF( KNUM.GT.JPLLMAX ) THEN 170 CALL INTLOG(JP_ERROR, 171 X 'WV2DINT: Number of latitudes in input lat/long grid = ',KNUM) 172 CALL INTLOG(JP_ERROR, 173 X 'WV2DINT: And is greater than allowed maximum = ',JPLLMAX) 174 WV2DINT = JPROUTINE + 2 175 GOTO 900 176 ENDIF 177C 178C Ensure working array dimensions are adequate for required output. 179C 180 IF( (KE_W.GT.JPNMOUT*2).OR.(KN_S.GT.JPNMOUT+1) ) THEN 181 CALL INTLOG(JP_ERROR, 182 X 'WV2DINT: Internal array dimensions are too small',JPQUIET) 183 CALL INTLOG(JP_ERROR, 184 X 'WV2DINT: for given lat/long output field.',JPQUIET) 185 WV2DINT = JPROUTINE + 3 186 GOTO 900 187 ENDIF 188C 189C Set up index for latitude lines in the input reduced lat/lon array 190C 191 INDEX(1) = 0 192 DO K = 2, KNUM 193 INDEX(K) = INDEX(K-1) + NUMPTS(K-1) 194 ENDDO 195C 196C Calculate latitudes and longitudes of output grid points. 197C 198C 199 DO K = 1, KN_S 200 LAT(K) = NORTH - FLOAT(K-1)*RESON 201 NUMNEW(K) = KE_W ! 0 202 ENDDO 203 204 DELAT = RNS/(KNUM-1) 205 206 DO k=1,KNUM 207 OLDLATS(K) = NORTH - FLOAT(K-1)*DELAT 208 END DO 209C 210 IF( ABS((KE_W*RESON)-360.0).LT.DEPS ) THEN 211 OLDEAST = WEST + KE_W*RESON 212 ELSE 213 OLDEAST = WEST + (KE_W-1)*RESON 214 ENDIF 215C 216C Initialise all points with 'missing data' indicator 217C 218 DO NEWROW = 1, KN_S 219 DO NEWCOL = 1, KE_W 220 NEWWAVE(NEWCOL,NEWROW) = PMISS 221 ENDDO 222 ENDDO 223C 224C --------------------------------------------------------------------- 225C* Section 2. Interpolation. 226C --------------------------------------------------------------------- 227C 228 200 CONTINUE 229C 230C Only calculate the indices on the first time through 231C 232 IF( (IFIRST.EQ.0 ) .OR. 233 X ( KNUM.NE.XKNUM ) .OR. 234 X ( KE_W.NE.XKE_W ) .OR. 235 X ( KN_S.NE.XKN_S ) .OR. 236 X ( RESON.NE.XRESON ).OR. 237 X ( NORTH.NE.XNORTH ).OR. 238 X ( WEST.NE.XWEST ) ) THEN 239 IRET = W251IDX(KNUM,NUMPTS,OLDLATS,WEST,OLDEAST, 240 X KN_S,NUMNEW,LAT,NEWIDX) 241 IFIRST = 1 242 XKNUM = KNUM 243 XKE_W = KE_W 244 XKN_S = KN_S 245 XRESON = RESON 246 XNORTH = NORTH 247 XWEST = WEST 248 ENDIF 249C 250C DELAT = 180.0/(KNUM-1) 251 DELAT = RNS/(KNUM-1) 252 NEXT = 0 253C 254 DO 220 NEWROW = 1, KN_S 255C 256C Set up the distance between new row and the old rows to N and S. 257C 258 NLAT = NINT((NORTH - LAT(NEWROW))/DELAT) + 1 259 I_N = NLAT 260 I_S = MIN(I_N+1,KNUM*2) 261 DIST_N = ((NORTH - DELAT*(I_N-1)) - LAT(NEWROW)) / DELAT 262 DIST_S = 1.0 - DIST_N 263C 264C Check if the new interpolated row lies between 2 old rows which 265C have data points 266C 267 IF( (NUMPTS(I_N).GT.0).AND.(NUMPTS(I_S).GT.0) ) THEN 268C 269C Yes, use the calculated indices 270C 271 DO 210 NEWCOL = 1, KE_W 272C 273 NEXT = NEXT + 1 274C 275 NEWWAVE(NEWCOL,NEWROW) = OLDWAVE(NEWIDX(NEXT)) 276C 277 210 CONTINUE 278 Else 279 NEXT = NEXT + KE_W 280 ENDIF 281C 282 220 CONTINUE 283C 284C --------------------------------------------------------------------- 285C* Section 9. Closedown. 286C --------------------------------------------------------------------- 287C 288 900 CONTINUE 289C 290 IF( WV2DINT.EQ.0 ) THEN 291 CALL INTLOG(JP_DEBUG, 292 X 'WV2DINT: Wave interpolation successful.',JPQUIET) 293 ELSE 294 CALL INTLOG(JP_DEBUG, 295 X 'WV2DINT: Wave interpolation failed.',JPQUIET) 296 ENDIF 297C 298 RETURN 299 END 300