MODULE W3SMCOMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Chris Bunney, UKMO | !/ | FORTRAN 90 | !/ | Last update : 27-Jun-2018 | !/ +-----------------------------------+ !/ !/ 18-Jan-2016 : Initial version. ( version 4.18 ) !/ (Chris Bunney, UK Met Office) !/ 28-Sep-2016 : Bug fix to EXO and EYO calcs for ( version 4.18 ) !/ whole domain output. (C. Bunney, UK Met Office) !/ 05-Oct-2016 : Bug fix to calculation of indicies ( version 4.18 ) !/ into regular grid for type 2 output. !/ (C. Bunney, UK Met Office) !/ 29-Sep-2017 : Revise calculation of indicies to ( version 4.18 ) !/ ensure selected cells fall inside !/ user selected areas. !/ (A. Saulter, UK Met Office) !/ 18-Apr-2018 : Added Type 3 and 4 SMC output: ( version 4.18 ) !/ re-gridding to an arbitrary regular !/ grid using pre-calculated indices and !/ weights (see WW3_SMCINT). Chris Bunney, UKMO !/ 20-Jun-2018 Updated routines to ensure that ( version 4.18 ) !/ directional fields are converted to !/ degrees (nautical convention) on output. !/ (Chris Bunney, UK Met Office) !/ 27-Jun-2018 Ported to V6 (Chris Bunney) ( version 6.05) !/ !/ Copyright 2009-2012 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! 1. Purpose : ! ! Service module for support of SMC => regular grid output ! ! 2. Method : ! ! For SMC grids, four types of output are possible: ! 1 - Flat grid (seapoint) output of SMC cells with associated ! cell size variables (cx and cy). Requires extra effort to ! plot as grid is not regular. ! ! 2 - Regular, uniformly gridded output to a specified output grid. ! This is achieved by area averaging of the SMC cells. The output ! grid will be aligned with the SMC grid cell edges which may ! result in the actual output grid being slightly different to ! the original user request. A land/sea mask is created by ! keeping track of the cell coverage and setting cells with ! total coverage of <50% UNDEF. ! ! 3 - Arbitrary regualar grid re-gridding using indices and weights ! generated from WW3_SMCINT. Uses the local gradient within the ! grid cell and distance between cell centres to interpolate. ! ! 4 - As type 3, but with no interpolation - value from surrounding ! SMC cell only. ! ! Note - directional fields are expected to be passed to into ! routines with units of RADIANS (CARTESIAN CONVENTION). ! They will be OUTPUT in units of DEGREES (NAUTICAL CONVENTION). ! USE W3GDATMD USE CONSTANTS USE W3ODATMD, ONLY: UNDEF PUBLIC ! Output grid definition REAL :: SXO, SYO, EXO, EYO ! First lat/lons ... REAL :: DXO, DYO ! ... grid size ... INTEGER :: NXO, NYO ! ... and number grid cells. ! Variables for SMC regridding (type 2 output): INTEGER :: SMCTYPE ! Type of SMC output INTEGER :: CELFAC ! Requested output cell scaling factor INTEGER, ALLOCATABLE :: XIDX(:), YIDX(:) ! Indices of SMC cells into regular grid INTEGER, ALLOCATABLE :: XSPAN(:), YSPAN(:) ! How many regualr grid cells SMC grid spans REAL, ALLOCATABLE :: WTS(:), COV(:,:) ! Regirdding weights and cell coverage (wet fraction) INTEGER, ALLOCATABLE :: MAPSMC(:,:) ! Regridded MAPSTA LOGICAL, ALLOCATABLE :: SMCMASK(:) ! Mask for type 1 output (flat array) INTEGER, ALLOCATABLE :: SMCIDX(:) ! Indices of SMC cells witin output grid domain INTEGER, ALLOCATABLE :: smccx(:), smccy(:) ! Cell size factors REAL :: dlat, dlon ! Base lat/lon sizes INTEGER :: cfac ! SMC scaling factor REAL :: NOVAL ! Fill value for seapoints with no value ! Variables for SMC nearest neighbour interpolation (type 3/4 output) INTEGER, ALLOCATABLE :: NNIDX(:,:) ! nearest neighbour smc point to regular grid REAL, ALLOCATABLE :: XDIST(:,:), YDIST(:,:) ! distance to nearest neighbour INTEGER :: NDSMC ! ww3_smcint file unit number ! Counters: INTEGER :: SMCNOUT, NSMC CONTAINS !-------------------------------------------------------------------------- SUBROUTINE SMC_INTERP() ! This subroutine generates index or mask values for extraction ! of SMC data to either a flat grid or regular lat/lon grid, ! depending on the type of SMC output grid selected: ! ! Type 1: Generates a mask for extracting only points from ! the user requested region. ! ! Type 2: Calculate interpolation indices and weights for ! regridding the irregular SMC grid onto a regular, ! uniformly spaced lat/lon grid. ! ! Chris Bunney, UK Met Office. 22-Oct-2015 IMPLICIT NONE ! Locals REAL :: CX0, CY0 ! SW corner of origin of grid REAL :: S0CHK INTEGER :: ISEA, mx, my, ixx, iyy, J REAL :: lat, lon J = 1 NSMC = 0 ! Determine smallest cell size factor: cfac = 2**(NRLv - 1) ! Get smallest SMC grid cells step size: dlat = SY / cfac dlon = SX / cfac ! SW Corner of grid origin cell: CX0 = X0 - SX / 2. CY0 = Y0 - SY / 2. ! Get start lat,lon (must be aligned with SMC grid edges). Use ! regular grid origins if SXO or SYO is -999.9 (use full grid): IF(ABS(SXO + 999.9) .LT. 1E-4) THEN SXO = CX0 ELSE S0CHK = CX0 + FLOOR((SXO - CX0) / SX) * SX ! Ensure first grid value falls within specified range ! NB the step is made at coarsest cell resolution IF (S0CHK .LT. SXO) THEN SXO = S0CHK + SX ELSE SXO = S0CHK ENDIF ENDIF IF(ABS(SYO + 999.9) .LT. 1E-4) THEN SYO = CY0 ELSE S0CHK = CY0 + FLOOR((SYO - CY0) / SY) * SY ! Ensure first grid value falls within specified range ! NB the step is made at coarsest cell resolution IF (S0CHK .LT. SYO) THEN SYO = S0CHK + SY ELSE SYO = S0CHK ENDIF ENDIF ! Use regular grid extents for last lat/lon if user ! specifies -999.9 for EXO/EYO (use full grid): IF(ABS(EXO + 999.9) .LT. 1E-4) THEN EXO = X0 + SX * (NX - 1) + 0.5 * SX !TRHC of last cell ! See below why the method below doesn't work: !EXO = CX0 + MAXVAL(IJKCel(1,:)) * dlon + dlon ENDIF IF(ABS(EYO + 999.9) .LT. 1E-4) THEN EYO = Y0 + SY * (NY-1) + 0.5 * SY !TRHC of last cell ! Method below doesn't work as the cell with the max val ! found might have a smaller value for CY, meaning another ! cell could extend further North. !EYO = CY0 + MAXVAL(IJKCel(2,:)) * dlat + dlat ENDIF ! Ouput grid cell dx/dy will be integer factor of smallest ! SMC grid cell size: DXO = dlon * celfac DYO = dlat * celfac ! Determine number of cells in grid: NXO = INT((EXO - SXO) / DXO) NYO = INT((EYO - SYO) / DYO) IF(SMCTYPE .EQ. 2) THEN ! Initialise all indices to "missing": XIDX(:) = -1 YIDX(:) = -1 ENDIF ! Loop over cell array and calculate regidding factors: DO ISEA=1, NSEA !/ARC ! For grids with Arctic region: make sure we don't double count !/ARC ! the overlapping boundary cells. Also, don't process the arctic !/ARC ! cell (which is always the last cell). !/ARC ! Note: NARC contains ALL the boundary cells (global + arctic). !/ARC ! whereas NGLO contains only the global boundary cells. !/ARC IF(ISEA .GT. NGLO-NBAC .AND. ISEA .LT. NSEA-NARC+1) CYCLE ! Get grid cell size: mx = IJKCel(3,ISEA) my = IJKCel(4,ISEA) ! Determine cell lat/lon (bottom left corner of cell) lon = CX0 + IJKCel(1,ISEA) * dlon lat = CY0 + IJKCel(2,ISEA) * dlat ! For output type 1 (seapoint array), just check whether ! cell centre is within specified domain range, and update ! output mask accordingly: IF( SMCTYPE .EQ. 1 ) THEN ! Ensure longitude ranges are aligned lon = lon + 0.5 * mx * dlon lat = lat + 0.5 * my * dlat IF(lon .LT. SXO) lon = lon + 360.0 IF(lon .GT. EXO) lon = lon - 360.0 ! Now check if it is within range of requested domain: IF(lon .GE. SXO .AND. lon .LE. EXO .AND. & lat .GE. SYO .AND. lat .LE. EYO ) THEN SMCMASK(ISEA) = .TRUE. SMCIDX(J) = ISEA J = J + 1 ENDIF CYCLE ENDIF ! SMCTYPE == 1 ! For output type 2 (area averaged regular grid), determine ! SMC grid cell location and coverage in output grid: ! Align lons IF(lon .LT. SXO) THEN lon = lon + 360. ELSE IF(lon .GT. EXO) THEN lon = lon - 360. ENDIF ! Find first SW cell in design grid: !!ixx = FLOOR((lon - SXO) / DXO) + 1 !!iyy = FLOOR((lat - SYO) / DYO) + 1 ! CB: Bugfix: In some grids, we end up with a precision issue using ! the above ixx/iyy calculations. Eg "(lon - SXO) / DXO" will end up ! being 4.9999992, rather than 5.0. Therefore, we add on 1/2 of the ! smallest SMC cell dlon/dlat values to ensure source grid cell ends ! up in the correct target grid cell: ixx = FLOOR((lon + 0.5*dlon - SXO) / DXO) + 1 iyy = FLOOR((lat + 0.5*dlat - SYO) / DYO) + 1 ! range check: IF(ixx .LE. 0 .OR. ixx .GT. NXO) THEN xidx(ISEA) = -1 yidx(ISEA) = -1 CYCLE ENDIF IF(iyy .LE. 0 .OR. iyy .GT. NYO) THEN xidx(ISEA) = -1 yidx(ISEA) = -1 CYCLE ENDIF XIDX(ISEA) = ixx YIDX(ISEA) = iyy NSMC = NSMC + 1 SMCIDX(NSMC) = ISEA ! find out how many cells it covers in the x/y directions: XSPAN(ISEA) = MAX(1, INT(mx / CELFAC)) YSPAN(ISEA) = MAX(1, INT(my / CELFAC)) ! Do a bit of error checking (non fatal - just produced warning): IF(XSPAN(ISEA) .GT. 1) THEN IF(sxo+(ixx-1)*dxo .NE. lon) THEN print*, 'Potential problem with grid cell span:' print*, xspan(ISEA), FLOAT(mx) / celfac print*, lon,lat print*, sxo+(ixx-1)*dxo,syo+iyy*dyo,dxo,dyo ENDIF ENDIF ! calc cell weight in relation to output grid: WTS(ISEA) = MIN(1., REAL(MIN(CELFAC, mx) * MIN(CELFAC, my)) / & (CELFAC**2)) ENDDO ! Reset SXO and SYO to be the cell-centre (currently cell SW edge): SXO = SXO + 0.5 * DXO SYO = SYO + 0.5 * DYO END SUBROUTINE SMC_INTERP !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE W3S2XY_SMCRG(S, XY) ! Regrid SMC data onto a regular grid using pre-calculated grid ! indices and weights. ! Chris Bunney, 02-Jul-2013 IMPLICIT NONE ! Input parameters: REAL, INTENT(IN) :: S(:) REAL, INTENT(OUT) :: XY(NXO,NYO) ! Local parameters INTEGER :: I, J, IX, IY, ISEA, ISMC ! Initialise coverage and output arrays: COV(:,:) = 0.0 XY(:,:) = 0.0 DO ISMC=1,NSMC ISEA = SMCIDX(ISMC) IF(S(ISEA) .EQ. UNDEF) CYCLE ! MDI ! Loop over number of spanned cells: DO I=0, XSPAN(ISEA) - 1 DO J=0, YSPAN(ISEA) - 1 IX = XIDX(ISEA) + I IY = YIDX(ISEA) + J ! Spans outside of grid? IF(IX .GT. NXO .OR. IY .GT. NYO) CYCLE ! Interpolate: XY(IX, IY) = XY(IX, IY) + S(ISEA) * WTS(ISEA) ! Keep track of how much of cell is (wet) covered: COV(IX, IY) = COV(IX, IY) + WTS(ISEA) ENDDO ENDDO ENDDO ! Create coastline by masking out areas with < 50% coverage: DO IX=1,NXO DO IY=1,NYO IF(MAPSMC(IX,IY) .EQ. 0) THEN ! Make land point XY(IX,IY) = UNDEF ELSE IF(COV(IX,IY) .LT. 0.5) THEN ! More than half of cell has UNDEF values - set to NOVAL: XY(IX,IY) = NOVAL ELSE IF(COV(IX,IY) .LT. 1.0) THEN ! If coverage < 1.0, scale values back to full cell coverage. ! Without this step, points around coast could end up with lower ! waveheights due to weights not summing to 1.0: XY(IX,IY) = XY(IX,IY) * ( 1.0 / COV(IX,IY) ) ENDIF ENDDO ENDDO RETURN END SUBROUTINE W3S2XY_SMCRG !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE W3S2XY_SMCRG_DIR(S, XY) ! Regrid SMC data from a directional field onto a regular grid using ! pre-calculated grid indices and weights. As W3S2XY_SMC, but ! decomposes the field into u/v components first to ensure proper ! area averaging of directional data. ! Chris Bunney, 02-Jul-2013 IMPLICIT NONE ! Input parameters: REAL, INTENT(IN) :: S(:) REAL, INTENT(OUT) :: XY(NXO,NYO) ! Local parameters INTEGER :: I, J, IX, IY, ISEA, ISMC REAL, ALLOCATABLE :: AUX1(:,:), AUX2(:,:) REAL :: COSS, SINS ! Initialise coverage and output arrays: ALLOCATE(AUX1(NXO,NYO),AUX2(NXO,NYO)) COV(:,:) = 0.0 XY(:,:) = 0.0 AUX1(:,:) = 0.0 AUX2(:,:) = 0.0 DO ISMC=1,NSMC ISEA = SMCIDX(ISMC) IF(S(ISEA) .EQ. UNDEF) CYCLE ! MDI COSS = COS(S(ISEA)) SINS = SIN(S(ISEA)) ! Loop over number of spanned cells: DO I=0, XSPAN(ISEA) - 1 DO J=0, YSPAN(ISEA) - 1 IX = XIDX(ISEA) + I IY = YIDX(ISEA) + J ! Spans outside of grid? IF(IX .GT. NXO .OR. IY .GT. NYO) CYCLE ! Interpolate: !XY(IX, IY) = XY(IX, IY) + S(ISEA) * WTS(ISEA) AUX1(IX, IY) = AUX1(IX, IY) + COSS * WTS(ISEA) AUX2(IX, IY) = AUX2(IX, IY) + SINS * WTS(ISEA) ! Keep track of how much of cell is (wet) covered: COV(IX, IY) = COV(IX, IY) + WTS(ISEA) ENDDO ENDDO ENDDO ! Create coastline by masking out areas with < 50% coverage: DO IX=1,NXO DO IY=1,NYO IF(MAPSMC(IX,IY) .EQ. 0) THEN ! Make land point XY(IX,IY) = UNDEF ELSE IF(COV(IX,IY) .LT. 0.5) THEN ! More than half of cell has UNDEF values - set to NOVAL XY(IX,IY) = NOVAL ELSE IF(COV(IX,IY) .LT. 1.0) THEN ! If coverage < 1.0, scale values back to full cell coverage. ! Without this step, points around coast could end up with lower ! waveheights due to weights not summing to 1.0: XY(IX,IY) = ATAN2(AUX2(IX,IY), AUX1(IX,IY)) XY(IX,IY) = MOD(630. - RADE * XY(IX,IY), 360. ) ELSE XY(IX,IY) = ATAN2(AUX2(IX,IY), AUX1(IX,IY)) XY(IX,IY) = MOD(630. - RADE * XY(IX,IY), 360. ) ENDIF ENDDO ENDDO RETURN END SUBROUTINE W3S2XY_SMCRG_DIR !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE MAPSTA_SMC() ! Calculates a new MATSTA using SMC grid cell averaging. ! Chris Bunney, 02-Jul-2013 IMPLICIT NONE ! Local parameters INTEGER :: I, J, IX, IY, IMX, IMY, ISEA ! Initialise coverage and output arrays: COV(:,:) = 0.0 MAPSMC(:,:) = 0 DO ISEA=1,NSEA IMX = MAPSF(ISEA,1) IMY = MAPSF(ISEA,2) IF(XIDX(ISEA) .EQ. -1) CYCLE ! Out of grid ! Loop over number of spanned cells: DO I=0, XSPAN(ISEA) - 1 DO J=0, YSPAN(ISEA) - 1 IX = XIDX(ISEA) + I IY = YIDX(ISEA) + J ! Spans outside of grid? IF(IX .GT. NXO .OR. IY .GT. NYO) CYCLE ! MAPSTA values: 0=Excluded, (+-)1=Sea, (+-2)=Input boundary ! We will just keep track of sea and non-sea points: IF(MAPSTA(IMY, IMX) .NE. 0) THEN ! Keep track of how much of cell is (wet) covered: COV(IX, IY) = COV(IX, IY) + WTS(ISEA) ENDIF ENDDO ENDDO ENDDO ! Create coastline by masking out areas with < 50% coverage: DO IX=1,NXO DO IY=1,NYO IF(COV(IX,IY) .LT. 0.5) THEN MAPSMC(IX, IY) = 0 ELSE MAPSMC(IX, IY) = 1 ENDIF ENDDO ENDDO RETURN END SUBROUTINE MAPSTA_SMC !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE READ_SMCINT() ! Reads the interpolation indices and distance weights from the ! smcint.ww3 file generated by ww3_smcint. !-------------------------------------------------------------------------- USE W3SERVMD, ONLY: EXTCDE IMPLICIT NONE ! Locals INTEGER :: IERR, I, J REAL :: PLATO, PLONO ! Not used yet....future version might allow ! output to a rotated pole grid... NDSMC = 50 OPEN(NDSMC, file='smcint.ww3', status='old', form='unformatted', iostat=ierr) IF(ierr .NE. 0) THEN WRITE(*,*) "ERROR! Failed to open smcint.ww3 for reading" CALL EXTCDE(1) ENDIF ! Header READ(NDSMC) NXO, NYO, SXO, SYO, DXO, DYO, PLONO, PLATO ALLOCATE(NNIDX(NXO,NYO), XDIST(NXO,NYO), YDIST(NXO,NYO)) ! Indices and weights: READ(NDSMC)((NNIDX(I,J), XDIST(I,J), YDIST(I,J),I=1,NXO),J=1,NYO) CLOSE(NDSMC) END SUBROUTINE READ_SMCINT !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE CALC_INTERP() ! Calculates the interpolation indices and weights for regirdding ! and SMC grid to an arbitrary regular grid. Index is that of ! SMC cell that contains output cell centre. Weights are the distance ! in metres between the output and SMC cell centres. ! A future version ~may~ allow for output grids to be on a rotated pole. !-------------------------------------------------------------------------- USE W3GDATMD, ONLY: CLATS USE CONSTANTS, ONLY : DERA, RADIUS !/RTD USE W3SERVMD, ONLY: W3LLTOEQ !/RTD USE W3GDATMD, ONLY: POLON, POLAT IMPLICIT NONE INTEGER :: IERR, I, J, ISEA, N, CFAC REAL :: mlon(NSEA), mlat(NSEA), olon(nxo,nyo), olat(nxo,nyo), & ang(nxo,nyo), lon, lat !/RTD REAL :: tmplon(nxo,nyo), tmplat(nxo,nyo) ! Determine smallest cell size factor: cfac = 2**(NRLv - 1) ! Get smallest SMC grid cells step size: dlat = SY / cfac dlon = SX / cfac ALLOCATE(xdist(nx,ny), ydist(ny,nx)) ! Model lat/lons: DO ISEA = 1,NSEA mlon(isea) = (X0-0.5*SX) + (IJKCel(1,ISEA) + 0.5 * IJKCel(3,ISEA)) * dlon mlat(isea) = (Y0-0.5*SY) + (IJKCel(2,ISEA) + 0.5 * IJKCel(4,ISEA)) * dlat ENDDO ! Generate output grid cell centres: DO I=1,NXO DO J=1,NYO olon(i,J) = SXO + (I-1) * DXO olat(i,J) = SYO + (J-1) * DYO ENDDO ENDDO !/RTD tmplat = olat !/RTD tmplon = olon !/RTD PRINT*,'Rotating coordinates' !/RTD CALL W3LLTOEQ ( tmplat, tmplon, olat, olon, & !/RTD ang, POLAT, POLON, NXO*NYO ) !/RTD PRINT*,'Rotating coordinates complete' ! Cycle over output grid points and find containing SMC cell: ! NOTE : BRUTE FORCE! NNIDX(:,:) = -1 DO I=1,NXO PRINT*,I,' of ',NXO DO J=1,NYO lon = olon(i,j) lat = olat(i,j) IF(lon .LT. X0 - SX / 2) lon = lon + 360.0 IF(lon .GT. (X0 + (NX-1) * SX) + 0.5 * SX) lon = lon - 360.0 DO ISEA=1,NSEA IF(mlon(ISEA) - 0.5 * IJKCel(3,ISEA) * dlon .LE. lon .AND. & mlon(ISEA) + 0.5 * IJKCel(3,ISEA) * dlon .GE. lon .AND. & mlat(ISEA) - 0.5 * IJKCel(4,ISEA) * dlat .LE. lat .AND. & mlat(ISEA) + 0.5 * IJKCel(4,ISEA) * dlat .GE. lat ) THEN ! Match! NNIDX(I,J) = ISEA xdist(I,J) = (lon - mlon(ISEA)) * DERA * RADIUS * CLats(ISEA) ydist(I,J) = (lat - mlat(ISEA)) * DERA * RADIUS EXIT ENDIF ENDDO ENDDO ENDDO END SUBROUTINE CALC_INTERP !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE W3S2XY_SMCNN(S, XY, DIRN) ! Fill regular grid using nearest SMC point data ! Directional fields (DIRN=True) will be assumed to be in radians ! and will be converted to degrees in nautical convention. !-------------------------------------------------------------------------- IMPLICIT NONE ! Input parameters: REAL, INTENT(IN) :: S(:) ! Inupt array REAL, INTENT(OUT) :: XY(NXO,NYO) ! Output data LOGICAL, INTENT(IN) :: DIRN ! Directional field? ! Local parameters INTEGER :: I, J, IX, IY, ISEA, ISMC DO IX = 1,NXO DO IY = 1,NYO ISEA = NNIDX(IX,IY) ! Nearest neighbour SMC point IF(ISEA .EQ. -1) THEN ! Land XY(IX,IY) = UNDEF ELSE IF(S(ISEA) .EQ. UNDEF) THEN ! Set undefined sea points to NOVAL XY(IX,IY) = NOVAL ELSE XY(IX,IY) = S(ISEA) IF(DIRN) THEN ! Convert direction fields to degrees nautical XY(IX,IY) = MOD(630. - RADE * XY(IX,IY), 360.0) ENDIF ENDIF ENDIF ENDDO ENDDO END SUBROUTINE W3S2XY_SMCNN !-------------------------------------------------------------------------- SUBROUTINE W3S2XY_SMCNN_INT(S, XY, DIRN) ! Fill regular grid using nearest SMC point data and interpolate ! output value based on local gradient and distance between grid ! cell centres. ! Directional fields (DIRN=True) will be assumed to be in radians ! and will be converted to degrees in nautical convention. !-------------------------------------------------------------------------- USE W3PSMCMD, ONLY: SMCGradn IMPLICIT NONE ! Input parameters: REAL, INTENT(IN) :: S(:) ! Input array REAL, INTENT(OUT) :: XY(NXO,NYO) ! Output array LOGICAL, INTENT(IN) :: DIRN ! Directional field? ! Locals INTEGER :: I, J, IX, IY, ISEA, ISMC REAL :: CVQ(-9:NSEA) REAL :: GrdX(NSEA), GrdY(NSEA) ! Calculate local gradients: CVQ(1:NSEA) = S ! Need to copy S into array with bounds starting at -9 CALL SMCGradn(CVQ, GrdX, GrdY, 0) ! Interpolate: DO IX = 1,NXO DO IY = 1,NYO ISEA = NNIDX(IX,IY) ! Nearest neighbour SMC point IF(ISEA .EQ. -1) THEN XY(IX,IY) = UNDEF ELSE ! Interpolate using local gradient and distance from cell centre: XY(IX,IY) = S(ISEA) + grdx(isea) * xdist(ix,iy) + grdy(isea) * ydist(ix,iy) IF(DIRN) THEN ! Convert direction fields to degrees nautical XY(IX,IY) = MOD(630. - RADE * XY(IX,IY), 360.0) ENDIF ENDIF ENDDO ENDDO END SUBROUTINE W3S2XY_SMCNN_INT !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- SUBROUTINE W3S2XY_SMC(S, XY, DIR) ! Entry point for SMC version of W3S2XY. ! Dispatches to regridding subroutine based on SMCTYPE. ! Optional DIR logical specifies whether field is a directional ! value; in which case it will be decomposed into u/v components ! prior to any interpolation. !-------------------------------------------------------------------------- IMPLICIT NONE REAL, INTENT(IN) :: S(:) REAL, INTENT(OUT) :: XY(NXO,NYO) LOGICAL, OPTIONAL :: DIR LOGICAL :: DIRN INTEGER :: ISEA IF(PRESENT(DIR)) THEN DIRN = DIR ELSE DIRN = .false. ENDIF IF(SMCTYPE .EQ. 1) THEN ! Flat sea point array XY(:,1) = PACK(S, SMCMASK) IF(DIRN) THEN ! Convert to nautical convention in degrees DO ISEA=1,NXO IF(XY(ISEA,1) .NE. UNDEF) THEN XY(ISEA,1) = MOD(630. - RADE * XY(ISEA,1), 360.) ENDIF ENDDO ENDIF ELSEIF(SMCTYPE .EQ. 2) THEN ! Regular gridded SMC cells IF(DIRN) THEN CALL W3S2XY_SMCRG_DIR(S, XY) ELSE CALL W3S2XY_SMCRG(S, XY) ENDIF ELSEIF(SMCTYPE .EQ. 3) THEN ! Regridded to arbitrary regular grid with interpolation CALL W3S2XY_SMCNN_INT(S, XY, DIRN) ELSEIF(SMCTYPE .EQ. 4) THEN ! Regridded to arbitrary regular grid - no interpolation CALL W3S2XY_SMCNN(S, XY, DIRN) ELSE WRITE(*,*) "Uknonwn SMC type!", SMCTYPE ! Unknown SMC type! STOP ENDIF END SUBROUTINE W3S2XY_SMC !-------------------------------------------------------------------------- END MODULE W3SMCOMD