w3smcomd.F90

Go to the documentation of this file.

 
!/
MODULE w3smcomd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |        Chris Bunney, UKMO         |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         21-Jul-2021 |
  !/                  +-----------------------------------+
  !/
  !/    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.
  !/
  USE w3gdatmd
  USE constants
  USE w3odatmd, ONLY: undef
 
  PUBLIC
 
  ! Output grid definition
  DOUBLE PRECISION     :: sxo
  DOUBLE PRECISION     :: syo
  DOUBLE PRECISION     :: exo
  DOUBLE PRECISION     :: eyo
  DOUBLE PRECISION     :: dxo
  DOUBLE PRECISION     :: dyo
  INTEGER              :: nxo
  INTEGER              :: nyo
 
  ! Variables for SMC regridding (type 2 output):
  INTEGER              :: smcotype
  INTEGER              :: celfac
  INTEGER, ALLOCATABLE :: xidx(:)
  INTEGER, ALLOCATABLE :: yidx(:)
  INTEGER, ALLOCATABLE :: xspan(:)
  INTEGER, ALLOCATABLE :: yspan(:)
  REAL, ALLOCATABLE    :: wts(:)
  REAL, ALLOCATABLE    :: cov(:,:)
  INTEGER, ALLOCATABLE :: mapsmc(:,:)
  LOGICAL, ALLOCATABLE :: smcmask(:)
  INTEGER, ALLOCATABLE :: smcidx(:)
 
  INTEGER, ALLOCATABLE :: smccx(:)
  INTEGER, ALLOCATABLE :: smccy(:)
  REAL                 :: dlat
  REAL                 :: dlon
  INTEGER              :: cfac
 
  REAL                 :: noval
 
  ! Variables for SMC nearest neighbour interpolation (type 3/4 output)
  INTEGER, ALLOCATABLE :: nnidx(:,:)
  REAL, ALLOCATABLE    :: xdist(:,:)
  REAL, ALLOCATABLE    :: ydist(:,:)
  INTEGER              :: ndsmc
 
  ! Counters:
  INTEGER              :: smcnout
  INTEGER              :: nsmc
 
CONTAINS
 
  !--------------------------------------------------------------------------
  SUBROUTINE smc_interp()
 
    IMPLICIT NONE
 
    ! Locals
    REAL    :: CX0, CY0   ! SW corner of origin of grid
    REAL    :: S0CHK, XSNAP, YSNAP
 
    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.
 
    ! Grid cell size to snap design grid to. Will be regular grid
    ! resolution for cellsize <= cfac, or cellsize for cellsize > cfac
    xsnap = sx
    ysnap = sy
    IF(celfac .gt. cfac) xsnap = celfac * dlon
    IF(celfac .gt. cfac) ysnap = celfac * dlat
 
    ! 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) / xsnap) * xsnap
      ! Ensure first grid value falls within specified range
      IF (s0chk .LT. sxo) THEN
        sxo = s0chk + xsnap
      ELSE
        sxo = s0chk
      ENDIF
    ENDIF
    IF(abs(syo + 999.9) .LT. 1e-4) THEN
      syo = cy0
    ELSE
      s0chk = cy0 + floor((syo - cy0) / ysnap) * ysnap
      ! Ensure first grid value falls within specified range
      IF (s0chk .LT. syo) THEN
        syo = s0chk + ysnap
      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 = cx0 + sx * nx ! TRHC of last cell
    ENDIF
    IF(abs(eyo + 999.9) .LT. 1e-4) THEN
      eyo = cy0 + sy * ny ! TRHC of last cell
    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 = nint((exo - sxo) / dxo)
    nyo = nint((eyo - syo) / dyo)
 
    IF(smcotype .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
      !          ! For grids with Arctic region: make sure we don't double count
      !          ! the overlapping boundary cells. Also, don't process the arctic
      !          ! cell (which is always the last cell).
      !          ! Note: NARC contains ALL the boundary cells (global + arctic).
      !          ! whereas NGLO contains only the global boundary cells.
      !          IF(ISEA .GT. NGLO-NBAC .AND. ISEA .LT. NSEA-NARC+1) CYCLE
      IF( arctc .AND. &
           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( smcotype .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 ! SMCOTYPE == 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.
      ENDIF
      IF(lon .GT. exo) THEN
        lon = lon - 360.
      ENDIF
 
      ! Find first SW cell in design grid:
      ! 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 (after
      ! integer trunction):
      ixx = floor((lon + 0.5*dlon - sxo) / dxo) + 1
      iyy = floor((lat + 0.5*dlat - syo) / dyo) + 1
 
      ! If we fall outside the left/bottom edge of the design grid,
      ! check for cases where the SMC cell has a lon or lat
      ! scaling factor > cfac (design grid is assumed to align
      ! its origin with cells of size cfac). For such cells,
      ! keep moving the left/bottom edge up by cfac until
      ! the SW corner (possibly) matches a design grid cell.
      IF(ixx .LE. 0 .AND. ixx + mx / celfac .GT. 0) THEN
        DO WHILE(mx .GT. cfac)
          mx = mx - cfac
          lon = lon + dlon * cfac
          ixx = floor((lon + 0.5*dlon - sxo) / dxo) + 1
          IF(ixx .GT. 0) EXIT ! Found cell lon-edge in design grid
        ENDDO
      ENDIF
      IF(iyy .LE. 0 .AND. iyy + my / celfac .GT. 0) THEN
        DO WHILE(my .GT. cfac)
          my = my - cfac
          lat = lat + dlat * cfac
          iyy = floor((lat + 0.5*dlat - syo) / dyo) + 1
          IF(iyy .GT. 0) EXIT ! Found cell lat-edge in design grid
        ENDDO
      ENDIF
 
      ! If SMC cell definitely out of design grid domain, then cycle.
      IF(ixx .LE. 0 .OR. ixx .GT. nxo .OR.                          &
           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(abs((sxo+(ixx-1)*dxo) - lon) .GT. dxo/100.0) THEN
          print*, 'Potential problem with SMC grid cell span:'
          print*, xspan(isea), float(mx) / celfac
          print*, lon,lat
          print*, sxo+(ixx-1)*dxo,syo+iyy*dyo,dxo,dyo
          print*, "diff:", (sxo+(ixx-1)*dxo) - lon
        ENDIF
      ENDIF
 
      ! calc cell weight in relation to output grid:
      wts(isea) = min(1., dble(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)
 
    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)
 
    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()
 
    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()
 
    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', convert=file_endian, 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()
 
    USE w3gdatmd, ONLY: clats
    USE constants, ONLY : dera, radius
#ifdef W3_RTD
    USE w3servmd, ONLY: w3lltoeq
    USE w3gdatmd, ONLY: polon, polat
#endif
 
    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
#ifdef W3_RTD
    REAL :: tmplon(nxo,nyo), tmplat(nxo,nyo)
#endif
 
    ! 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
 
#ifdef W3_RTD
    tmplat = olat
    tmplon = olon
    print*,'Rotating coordinates'
    CALL w3lltoeq ( tmplat, tmplon, olat, olon,                       &
         ang, polat, polon, nxo*nyo )
    print*,'Rotating coordinates complete'
#endif
 
    ! 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)
 
    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)
 
    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)
 
    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(smcotype .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(smcotype .EQ. 2) THEN
      ! Regular gridded SMC cells
      IF(dirn) THEN
        CALL w3s2xy_smcrg_dir(s, xy)
      ELSE
        CALL w3s2xy_smcrg(s, xy)
      ENDIF
    ELSEIF(smcotype .EQ. 3) THEN
      ! Regridded to arbitrary regular grid with interpolation
      CALL w3s2xy_smcnn_int(s, xy, dirn)
    ELSEIF(smcotype .EQ. 4) THEN
      ! Regridded to arbitrary regular grid - no interpolation
      CALL w3s2xy_smcnn(s, xy, dirn)
    ELSE
      WRITE(*,*) "Uknonwn SMC type!", smcotype
      ! Unknown SMC type!
      stop
    ENDIF
 
  END SUBROUTINE w3s2xy_smc
  !--------------------------------------------------------------------------
 
END MODULE w3smcomd