WAVEWATCH III  beta 0.0.1
Functions/Subroutines
wmscrpmd Module Reference

Routines to determine and process grid dependencies in the multi-grid wave model. More...

Functions/Subroutines

subroutine scrip_wrapper (ID_SRC, ID_DST, MAPSTA_SRC, MAPST2_SRC, FLAGLL, GRIDSHIFT, L_MASTER, L_READ, L_TEST)
 Compute grid information required by SCRIP. More...
 
subroutine get_scrip_info_unstructured (ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
 Compute grid arrays for scrip for a specific unstructured grid. More...
 
subroutine get_scrip_info_structured (ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
 Compute grid arrays needed for scrip for a specific structured grid. More...
 
subroutine get_scrip_info (ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
 Compute grid for scrip for a specific structured grid. More...
 
subroutine scrip_info_renormalization (GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK, CONV_DX, CONV_DY, OFFSET, GRIDSHIFT)
 Rescale according to whether the grid is spherical or not. More...
 
subroutine triang_indexes (I, INEXT, IPREV)
 Desc not available. More...
 
subroutine get_unstructured_vertex_degree (MNP, MNE, TRIGP, TRIGINCD)
 This function returns the list of incidences. More...
 
subroutine get_boundary (MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV, NEIGHBOR_NEXT)
 Returns neighbor of a boundary node. More...
 
subroutine fix_periodcity (PT)
 Adjust element longitude coordinates for elements straddling the dateline with distance of ~360 degrees. More...
 

Detailed Description

Routines to determine and process grid dependencies in the multi-grid wave model.

Author
E. Rogers
M. Dutour
A. Roland
F. Ardhuin
Date
10-Dec-2014

Function/Subroutine Documentation

◆ fix_periodcity()

subroutine wmscrpmd::fix_periodcity ( real*8, dimension(3,2), intent(inout)  PT)

Adjust element longitude coordinates for elements straddling the dateline with distance of ~360 degrees.

Detect if element has nodes on both sides of dateline and adjust coordinates so that all nodes have the same sign.

Parameters
[in,out]PT
Author
Steven Brus
Ali Abdolali
Date
21-May-2020

Definition at line 1711 of file wmscrpmd.F90.

1711 
1712  !/
1713  !/ +-----------------------------------+
1714  !/ | WAVEWATCH III NOAA/NCEP |
1715  !/ | Steven Brus |
1716  !/ | Ali Abdolali |
1717  !/ | FORTRAN 90 |
1718  !/ | Last update : 21-May-2020 |
1719  !/ +-----------------------------------+
1720  !/
1721  !/ 21-May-2020 : Origination. ( version 6.07 )
1722  !/
1723  !/
1724  ! 1. Purpose :
1725  !
1726  ! Adjust element longitude coordinates for elements straddling the
1727  ! dateline with distance of ~360 degrees
1728  !
1729  ! 2. Method :
1730  !
1731  ! Detect if element has nodes on both sides of dateline and adjust
1732  ! coordinates so that all nodes have the same sign
1733  !
1734  ! 3. Parameters :
1735  !
1736  ! Parameter list
1737  ! ----------------------------------------------------------------
1738  IMPLICIT NONE
1739  real*8, INTENT(INOUT) :: pt(3,2)
1740  ! ----------------------------------------------------------------
1741  !
1742  ! Local variables.
1743  ! ----------------------------------------------------------------
1744  INTEGER :: I
1745  INTEGER :: R1GT180, R2GT180, R3GT180
1746  ! ----------------------------------------------------------------
1747  !
1748  ! 4. Subroutines used :
1749  !
1750 
1751  ! 5. Called by :
1752  !
1753  ! Name Type Module Description
1754  ! ----------------------------------------------------------------
1755  ! GET_SCRIP_INFO_UNSTRUCTURED Subr. WMSCRPMD Element center calculation
1756  ! GET_SCRIP_INFO Subr. WMSCRPMD Check signs
1757  ! ----------------------------------------------------------------
1758  !
1759  ! 6. Error messages :
1760  !
1761  ! None.
1762  !
1763  ! 7. Remarks :
1764  !
1765  ! 8. Structure :
1766  !
1767  ! 9. Switches :
1768  !
1769  ! 10. Source code :
1770  !/ ------------------------------------------------------------------- /
1771 
1772  r1gt180 = merge(1, 0, abs(pt(3,1)-pt(2,1)).GT.180)
1773  r2gt180 = merge(1, 0, abs(pt(1,1)-pt(3,1)).GT.180)
1774  r3gt180 = merge(1, 0, abs(pt(2,1)-pt(1,1)).GT.180)
1775  ! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the
1776  ! dateline
1777  ! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole
1778  ! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline
1779 
1780  IF ( r1gt180 + r2gt180 == 2 ) THEN
1781  pt(3,1)=pt(3,1)-sign(360.0d0,(pt(3,1)-pt(2,1)))
1782  ELSE IF ( r2gt180 + r3gt180 == 2 ) THEN
1783  pt(1,1)=pt(1,1)-sign(360.0d0,(pt(1,1)-pt(2,1)))
1784  ELSE IF ( r1gt180 + r3gt180 == 2 ) THEN
1785  pt(2,1)=pt(2,1)-sign(360.0d0,(pt(2,1)-pt(3,1)))
1786  ENDIF
1787 
1788  RETURN

Referenced by get_scrip_info(), and get_scrip_info_unstructured().

◆ get_boundary()

subroutine wmscrpmd::get_boundary ( integer, intent(in)  MNP,
integer, intent(in)  MNE,
integer, dimension(3,mne), intent(in)  TRIGP,
integer, dimension(mnp), intent(inout)  IOBP,
integer, dimension(mnp), intent(inout)  NEIGHBOR_PREV,
integer, dimension(mnp), intent(inout)  NEIGHBOR_NEXT 
)

Returns neighbor of a boundary node.

If a node belong to a boundary, the function returns the neighbor of this point on one side. If the point is interior then the value 0 is set.

Parameters
[in]MNPNumber of nodes.
[in]MNENumber of triangles.
[in]TRIGPList of nodes.
[in,out]IOBP
[in,out]NEIGHBOR_PREV
[in,out]NEIGHBOR_NEXT
Author
M. Dutour
A. Roland
Date
10-Dec-2014

Definition at line 1502 of file wmscrpmd.F90.

1502  !/ +-----------------------------------+
1503  !/ | WAVEWATCH III |
1504  !/ | M. Dutour, A. Roland |
1505  !/ | FORTRAN 90 |
1506  !/ | Last update : 10-Dec-2014 !
1507  !/ +-----------------------------------+
1508  !/
1509  ! Written:
1510  !
1511  ! 20-Feb-2012
1512  !/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
1513  !
1514  ! Author:
1515  !
1516  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
1517  !
1518  ! Parameters:
1519  ! Input:
1520  ! MNP: number of nodes
1521  ! TRIGP: list of nodes
1522  ! MNE: number of triangles
1523  ! Output:
1524  ! NEIGHBOR
1525  !
1526  ! Description:
1527  ! if a node belong to a boundary, the function
1528  ! returns the neighbor of this point on one side.
1529  ! if the point is interior then the value 0 is set.
1530  !
1531  USE w3servmd, ONLY: extcde
1532  IMPLICIT NONE
1533 
1534  INTEGER, INTENT(IN) :: MNP, MNE, TRIGP(3,MNE)
1535  INTEGER, INTENT(INOUT) :: IOBP(MNP)
1536  INTEGER, INTENT(INOUT) :: NEIGHBOR_PREV(MNP)
1537  INTEGER, INTENT(INOUT) :: NEIGHBOR_NEXT(MNP)
1538 
1539  INTEGER, POINTER :: STATUS(:)
1540  INTEGER, POINTER :: COLLECTED(:)
1541  INTEGER, POINTER :: NEXTVERT(:)
1542  INTEGER, POINTER :: PREVVERT(:)
1543 
1544  INTEGER :: IE, I, IP, IP2, IP3
1545  INTEGER :: ISFINISHED, INEXT, IPREV
1546  INTEGER :: IPNEXT, IPPREV, ZNEXT, ZPREV
1547 
1548  ALLOCATE(status(mnp), stat=istat)
1549  check_alloc_status( istat )
1550  ALLOCATE(collected(mnp), stat=istat)
1551  check_alloc_status( istat )
1552  ALLOCATE(prevvert(mnp), stat=istat)
1553  check_alloc_status( istat )
1554  ALLOCATE(nextvert(mnp), stat=istat)
1555  check_alloc_status( istat )
1556  iobp = 0
1557  neighbor_next = 0
1558  neighbor_prev = 0
1559 
1560  ! Now computing the next items
1561  status = 0
1562  nextvert = 0
1563  prevvert = 0
1564  DO ie=1,mne
1565  DO i=1,3
1566  CALL triang_indexes(i, inext, iprev)
1567  ip=trigp(i,ie)
1568  ipnext=trigp(inext,ie)
1569  ipprev=trigp(iprev,ie)
1570  IF (status(ip).EQ.0) THEN
1571  status(ip)=1
1572  prevvert(ip)=ipprev
1573  nextvert(ip)=ipnext
1574  END IF
1575  END DO
1576  END DO
1577  status(:)=0
1578  DO
1579  collected(:)=0
1580  DO ie=1,mne
1581  DO i=1,3
1582  CALL triang_indexes(i, inext, iprev)
1583  ip=trigp(i,ie)
1584  ipnext=trigp(inext,ie)
1585  ipprev=trigp(iprev,ie)
1586  IF (status(ip).EQ.0) THEN
1587  znext=nextvert(ip)
1588  IF (znext.EQ.ipprev) THEN
1589  collected(ip)=1
1590  nextvert(ip)=ipnext
1591  IF (nextvert(ip).EQ.prevvert(ip)) THEN
1592  status(ip)=1
1593  END IF
1594  END IF
1595  END IF
1596  END DO
1597  END DO
1598 
1599  isfinished=1
1600  DO ip=1,mnp
1601  IF ((collected(ip).EQ.0).AND.(status(ip).EQ.0)) THEN
1602  status(ip)=-1
1603  neighbor_next(ip)=nextvert(ip)
1604  END IF
1605  IF (status(ip).EQ.0) THEN
1606  isfinished=0
1607  END IF
1608  END DO
1609  IF (isfinished.EQ.1) THEN
1610  EXIT
1611  END IF
1612  END DO
1613 
1614  ! Now computing the prev items
1615  status = 0
1616  nextvert = 0
1617  prevvert = 0
1618  DO ie=1,mne
1619  DO i=1,3
1620  CALL triang_indexes(i, inext, iprev)
1621  ip=trigp(i,ie)
1622  ipnext=trigp(inext,ie)
1623  ipprev=trigp(iprev,ie)
1624  IF (status(ip).EQ.0) THEN
1625  status(ip)=1
1626  prevvert(ip)=ipprev
1627  nextvert(ip)=ipnext
1628  END IF
1629  END DO
1630  END DO
1631  status(:)=0
1632  DO
1633  collected(:)=0
1634  DO ie=1,mne
1635  DO i=1,3
1636  CALL triang_indexes(i, inext, iprev)
1637  ip=trigp(i,ie)
1638  ipnext=trigp(inext,ie)
1639  ipprev=trigp(iprev,ie)
1640  IF (status(ip).EQ.0) THEN
1641  zprev=prevvert(ip)
1642  IF (zprev.EQ.ipnext) THEN
1643  collected(ip)=1
1644  prevvert(ip)=ipprev
1645  IF (prevvert(ip).EQ.nextvert(ip)) THEN
1646  status(ip)=1
1647  END IF
1648  END IF
1649  END IF
1650  END DO
1651  END DO
1652 
1653  isfinished=1
1654  DO ip=1,mnp
1655  IF ((collected(ip).EQ.0).AND.(status(ip).EQ.0)) THEN
1656  status(ip)=-1
1657  neighbor_prev(ip)=prevvert(ip) ! new code
1658  END IF
1659  IF (status(ip).EQ.0) THEN
1660  isfinished=0
1661  END IF
1662  END DO
1663  IF (isfinished.EQ.1) THEN
1664  EXIT
1665  END IF
1666  END DO
1667  ! Now making checks
1668  DO ip=1,mnp
1669  ip2=neighbor_next(ip)
1670  IF (ip2.GT.0) THEN
1671  ip3=neighbor_prev(ip2)
1672  IF (abs(ip3 - ip).GT.0) THEN
1673  WRITE(*,*) 'IP=', ip, ' IP2=', ip2, ' IP3=', ip3
1674  WRITE(*,*) 'We have a dramatic inconsistency'
1675  stop 'wmscrpmd, case 3'
1676  END IF
1677  END IF
1678  END DO
1679  ! Now assigning the boundary IOBP array
1680  DO ip=1,mnp
1681  IF (status(ip).EQ.-1 .AND. iobp(ip) .EQ. 0) THEN
1682  iobp(ip)=1
1683  END IF
1684  END DO
1685 
1686  DEALLOCATE(status, stat=istat)
1687  check_dealloc_status( istat )
1688  DEALLOCATE(collected, stat=istat)
1689  check_dealloc_status( istat )
1690  DEALLOCATE(nextvert, stat=istat)
1691  check_dealloc_status( istat )
1692  DEALLOCATE(prevvert, stat=istat)
1693  check_dealloc_status( istat )
1694 

References w3servmd::extcde(), and triang_indexes().

Referenced by get_scrip_info_unstructured().

◆ get_scrip_info()

subroutine wmscrpmd::get_scrip_info ( integer, intent(in)  ID_GRD,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LON,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LAT,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LON,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LAT,
logical, dimension(:), intent(out), allocatable  GRID_MASK,
integer, dimension(:), intent(out), allocatable  GRID_DIMS,
integer, intent(out)  GRID_SIZE,
integer, intent(out)  GRID_CORNERS,
integer, intent(out)  GRID_RANK 
)

Compute grid for scrip for a specific structured grid.

This is adapted from Erick Rogers code by making it cleaner.

Parameters
[in]ID_GRD
[out]GRID_CENTER_LON
[out]GRID_CENTER_LAT
[out]GRID_CORNER_LON
[out]GRID_CORNER_LAT
[out]GRID_MASK
[out]GRID_DIMS
[out]GRID_SIZE
[out]GRID_CORNERS
[out]GRID_RANK
Author
M. Dutour
A. Roland
Date
20-Feb-2012

Definition at line 1165 of file wmscrpmd.F90.

1165  ! 1. Original author :
1166  !
1167  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
1168  !
1169  ! 2. Last update :
1170  !
1171  ! See revisions.
1172  !
1173  ! 3. Revisions :
1174  !
1175  ! 20-Feb-2012 : Origination, this is adapted from Erick Rogers
1176  ! code by splitting the code into sections.
1177  !
1178  ! 4. Copyright :
1179  !
1180  ! 5. Purpose :
1181  !
1182  ! Compute grid for scrip for a specific structured grid.
1183  ! This is adapted from Erick Rogers code by making it cleaner.
1184  !
1185  ! 6. Method :
1186  !
1187  ! 7. Parameters, Variables and types :
1188  !
1189  ! 8. Called by :
1190  !
1191  ! Subroutine scrip_wrapper
1192  !
1193  ! 9. Subroutines and functions used :
1194  !
1195  ! 10. Error messages:
1196  !
1197  ! 11. Remarks :
1198  !
1199  ! 12. Structure :
1200  !
1201  ! 13. Switches :
1202  !
1203  ! 14. Source code :
1204  USE w3servmd, ONLY: extcde
1205  USE w3gdatmd, ONLY : grids, ungtype
1206  IMPLICIT NONE
1207  INTEGER, INTENT(IN) :: ID_GRD
1208  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
1209  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
1210  LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
1211  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
1212  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
1213  INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
1214  INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
1215  real*8 :: dlon1, dlat1, dlon2, dlat2, thedet
1216  INTEGER :: I, J
1217  INTEGER :: IC, JC, IP, CHECKSIGNS, NBPLUS, NBMINUS, NBZERO
1218  INTEGER :: PRINTDATA, PRINTMINMAX
1219  real*8 :: minlon, minlat, maxlon, maxlat
1220  real*8 :: minloncorner, maxloncorner, minlatcorner, maxlatcorner
1221  real*8 :: pt(3,2)
1222  IF (grids(id_grd)%GTYPE .EQ. ungtype) THEN
1223  CALL get_scrip_info_unstructured (id_grd, &
1224  & grid_center_lon, grid_center_lat, &
1225  & grid_corner_lon, grid_corner_lat, grid_mask, &
1226  & grid_dims, grid_size, grid_corners, grid_rank)
1227  ELSE
1228  CALL get_scrip_info_structured (id_grd, &
1229  & grid_center_lon, grid_center_lat, &
1230  & grid_corner_lon, grid_corner_lat, grid_mask, &
1231  & grid_dims, grid_size, grid_corners, grid_rank)
1232  END IF
1233  checksigns=1
1234  IF (checksigns.EQ.1) THEN
1235  nbplus=0
1236  nbminus=0
1237  nbzero=0
1238  DO ip=1,grid_size
1239  DO ic=1,grid_corners
1240  IF (ic.EQ.grid_corners) THEN
1241  jc=1
1242  ELSE
1243  jc=ic+1
1244  END IF
1245 
1246  pt(1,1) = grid_center_lon(ip)
1247  pt(1,2) = grid_center_lat(ip)
1248  pt(2,1) = grid_corner_lon(ic,ip)
1249  pt(2,2) = grid_corner_lat(ic,ip)
1250  pt(3,1) = grid_corner_lon(jc,ip)
1251  pt(3,2) = grid_corner_lat(jc,ip)
1252 
1253  CALL fix_periodcity(pt)
1254 
1255  dlon1=pt(2,1)-pt(1,1)
1256  dlon2=pt(3,1)-pt(1,1)
1257  dlat1=pt(2,2)-pt(1,2)
1258  dlat2=pt(3,2)-pt(1,2)
1259 
1260  thedet=dlon1*dlat2 - dlon2*dlat1
1261  IF (thedet.GT.1d-8) THEN
1262  nbplus=nbplus+1
1263  ELSE IF (thedet.LT.-1d-8) THEN
1264  nbminus=nbminus+1
1265  ELSE
1266  nbzero=nbzero+1
1267  END IF
1268  END DO
1269  END DO
1270 
1271  WRITE(*,*) 'SI nbPlus=', nbplus, ' nbMinus=', nbminus, ' nbZero=', nbzero
1272 
1273  END IF

References w3servmd::extcde(), fix_periodcity(), get_scrip_info_structured(), get_scrip_info_unstructured(), w3gdatmd::grids, and w3gdatmd::ungtype.

Referenced by scrip_wrapper().

◆ get_scrip_info_structured()

subroutine wmscrpmd::get_scrip_info_structured ( integer, intent(in)  ID_GRD,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LON,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LAT,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LON,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LAT,
logical, dimension(:), intent(out), allocatable  GRID_MASK,
integer, dimension(:), intent(out), allocatable  GRID_DIMS,
integer, intent(out)  GRID_SIZE,
integer, intent(out)  GRID_CORNERS,
integer, intent(out)  GRID_RANK 
)

Compute grid arrays needed for scrip for a specific structured grid.

This is adapted from Erick Rogers original code by splitting the original scrip_wrapper function.

Parameters
[in]ID_GRD
[out]GRID_CENTER_LON
[out]GRID_CENTER_LAT
[out]GRID_CORNER_LON
[out]GRID_CORNER_LAT
[out]GRID_MASK
[out]GRID_DIMS
[out]GRID_SIZE
[out]GRID_CORNERS
[out]GRID_RANK
Author
M. Dutour
A. Roland
Date
10-Dec-2014

Definition at line 974 of file wmscrpmd.F90.

974  !/ +-----------------------------------+
975  !/ | WAVEWATCH III |
976  !/ | M. Dutour, A. Roland |
977  !/ | FORTRAN 90 |
978  !/ | Last update : 10-Dec-2014 !
979  !/ +-----------------------------------+
980  !/
981  ! 1. Original author :
982  !
983  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
984  !
985  ! 2. Last update :
986  !
987  ! See revisions.
988  !
989  ! 3. Revisions :
990  !
991  ! 20-Feb-2012 : Origination, this is adapted from Erick Rogers
992  ! code by splitting the code into sections.
993  !/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
994  !
995  ! 4. Copyright :
996  !
997  ! 5. Purpose :
998  !
999  ! Compute grid arrays needed for scrip for a specific structured grid.
1000  ! This is adapted from Erick Rogers original code by splitting
1001  ! the original scrip_wrapper function
1002  !
1003  ! 6. Method :
1004  !
1005  ! 7. Parameters, Variables and types :
1006  !
1007  ! 8. Called by :
1008  !
1009  ! Subroutine get_scrip_info
1010  !
1011  ! 9. Subroutines and functions used :
1012  !
1013  ! 10. Error messages:
1014  !
1015  ! 11. Remarks :
1016  !
1017  ! 12. Structure :
1018  !
1019  ! 13. Switches :
1020  !
1021  ! 14. Source code :
1022  USE w3servmd, ONLY: extcde
1023  USE w3gdatmd, ONLY : grids
1024  USE scrip_constants, ONLY : half
1025  IMPLICIT NONE
1026  INTEGER, INTENT(IN) :: ID_GRD
1027  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
1028  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
1029  LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
1030  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
1031  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
1032  INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
1033  INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
1034  !
1035  real*8, ALLOCATABLE :: xin_grd(:,:), yin_grd(:,:)
1036  real*8, ALLOCATABLE :: dxdp_grd(:,:), dxdq_grd(:,:)
1037  real*8, ALLOCATABLE :: dydp_grd(:,:), dydq_grd(:,:)
1038  INTEGER :: N1, N2, NI, NJ
1039  INTEGER :: IREC, J, I
1040  grid_rank=2
1041  n1=SIZE(grids(id_grd)%XGRD,1)
1042  n2=SIZE(grids(id_grd)%XGRD,2)
1043  ALLOCATE(xin_grd(n1,n2), stat=istat)
1044  check_alloc_status( istat )
1045  ALLOCATE(yin_grd(n1,n2), stat=istat)
1046  check_alloc_status( istat )
1047  ALLOCATE(dxdp_grd(n1,n2), stat=istat)
1048  check_alloc_status( istat )
1049  ALLOCATE(dxdq_grd(n1,n2), stat=istat)
1050  check_alloc_status( istat )
1051  ALLOCATE(dydp_grd(n1,n2), stat=istat)
1052  check_alloc_status( istat )
1053  ALLOCATE(dydq_grd(n1,n2), stat=istat)
1054  check_alloc_status( istat )
1055 
1056  xin_grd=dble(grids(id_grd)%XGRD)
1057  yin_grd=dble(grids(id_grd)%YGRD)
1058  dxdp_grd=dble(grids(id_grd)%DXDP)
1059  dxdq_grd=dble(grids(id_grd)%DXDQ)
1060  dydp_grd=dble(grids(id_grd)%DYDP)
1061  dydq_grd=dble(grids(id_grd)%DYDQ)
1062 
1063  ALLOCATE(grid_dims(grid_rank), stat=istat)
1064  check_alloc_status( istat )
1065  grid_dims(1)=n2
1066  ni=n2
1067  grid_dims(2)=n1
1068  nj=n1
1069 
1070  grid_size=ni*nj ! hardwired: logically rectangular grid
1071  grid_corners=4 ! hardwired: each cell has 4 corners
1072 
1073  !.....notes: unfortunately, scrip only works for spherical coordinates.
1074  ! thus, if we want to have a multi-grid case in meters, we have to
1075  ! fake it. fortunately, it should be pretty rare to have a multi-grid
1076  ! case in meters.
1077 
1078  ALLOCATE(grid_center_lon(ni*nj), stat=istat)
1079  check_alloc_status( istat )
1080  ALLOCATE(grid_center_lat(ni*nj), stat=istat)
1081  check_alloc_status( istat )
1082  ALLOCATE(grid_corner_lon(4,ni*nj), stat=istat)
1083  check_alloc_status( istat )
1084  ALLOCATE(grid_corner_lat(4,ni*nj), stat=istat)
1085  check_alloc_status( istat )
1086  ALLOCATE(grid_mask(ni*nj), stat=istat)
1087  check_alloc_status( istat )
1088 
1089  !.....notes: this "gridshift" variable is included because SCRIP sometimes
1090  ! has trouble when grids cell locations are identical between
1091  ! the two grids. Thus we apply this to one of the two grids.
1092 
1093  !.....notes: The following block of code could be converted to a subroutine.
1094  ! Since it is called twice, this would save a little space.
1095  irec=0
1096  DO j=1,nj
1097  DO i=1,ni
1098  irec=irec+1
1099  grid_center_lon(irec)=xin_grd(j,i)
1100  grid_center_lat(irec)=yin_grd(j,i)
1101  grid_mask(irec)=.true.
1102 
1103  !..notes: normally, we'd apply the mask like this:
1104  ! if(abs(mapsta_src(j,i)).eq.1)then
1105  ! grid1_mask(irec)=.true.
1106  ! else
1107  ! grid1_mask(irec)=.false.
1108  ! endif
1109  !..but unfortunately, WMGHGH needs information about the overlaying high-res
1110  ! cells, even those that are masked, for calculating
1111  ! NRL, NR0, NR1, NR2.
1112 
1113  !...........corner 1 : halfway to i-1,j-1
1114  grid_corner_lon(1,irec)=grid_center_lon(irec)- &
1115  & half*dxdp_grd(j,i)-half*dxdq_grd(j,i)
1116  grid_corner_lat(1,irec)=grid_center_lat(irec)- &
1117  & half*dydp_grd(j,i)-half*dydq_grd(j,i)
1118 
1119  !...........corner 2: halfway to i+1,j-1
1120  grid_corner_lon(2,irec)=grid_center_lon(irec)+ &
1121  & half*dxdp_grd(j,i)-half*dxdq_grd(j,i)
1122  grid_corner_lat(2,irec)=grid_center_lat(irec)+ &
1123  & half*dydp_grd(j,i)-half*dydq_grd(j,i)
1124 
1125  !...........corner 3: halfway to i+1,j+1
1126  grid_corner_lon(3,irec)=grid_center_lon(irec)+ &
1127  & half*dxdp_grd(j,i)+half*dxdq_grd(j,i)
1128  grid_corner_lat(3,irec)=grid_center_lat(irec)+ &
1129  & half*dydp_grd(j,i)+half*dydq_grd(j,i)
1130 
1131  !...........corner 4: halfway to i-1,j+1
1132  grid_corner_lon(4,irec)=grid_center_lon(irec)- &
1133  & half*dxdp_grd(j,i)+half*dxdq_grd(j,i)
1134  grid_corner_lat(4,irec)=grid_center_lat(irec)- &
1135  & half*dydp_grd(j,i)+half*dydq_grd(j,i)
1136  END DO
1137  END DO

References w3servmd::extcde(), w3gdatmd::grids, and scrip_constants::half.

Referenced by get_scrip_info().

◆ get_scrip_info_unstructured()

subroutine wmscrpmd::get_scrip_info_unstructured ( integer, intent(in)  ID_GRD,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LON,
real*8, dimension(:), intent(out), allocatable  GRID_CENTER_LAT,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LON,
real*8, dimension(:,:), intent(out), allocatable  GRID_CORNER_LAT,
logical, dimension(:), intent(out), allocatable  GRID_MASK,
integer, dimension(:), intent(out), allocatable  GRID_DIMS,
integer, intent(out)  GRID_SIZE,
integer, intent(out)  GRID_CORNERS,
integer, intent(out)  GRID_RANK 
)

Compute grid arrays for scrip for a specific unstructured grid.

For interior vertices, we select for every triangle the barycenter of the triangle. So to every node contained in N triangles we associate a domain with N corners. Every one of those corners is contained in 3 different domains.

For nodes that are on the boundary, we have to proceed differently. For every such node, we have NEIGHBOR_PREV and NEIGHBOR_NEXT which are the neighbor on each side of the boundary. We put a corner on the middle of the edge. We also put a corner on the node itself.

Note that instead of taking barycenters, we could have taken Voronoi vertices, but this carries danger since Voronoi vertices can be outside of the triangle. And it leaves open how to treat the boundary.

Parameters
[in]ID_GRD
[out]GRID_CENTER_LON
[out]GRID_CENTER_LAT
[out]GRID_CORNER_LON
[out]GRID_CORNER_LAT
[out]GRID_MASK
[out]GRID_DIMS
[out]GRID_SIZE
[out]GRID_CORNERS
[out]GRID_RANK
Author
M. Dutour
A. Roland
Date
10-Dec-2014

Definition at line 551 of file wmscrpmd.F90.

551  !/ +-----------------------------------+
552  !/ | WAVEWATCH III |
553  !/ | M. Dutour, A. Roland |
554  !/ | FORTRAN 90 |
555  !/ | Last update : 10-Dec-2014 !
556  !/ +-----------------------------------+
557  !/
558  ! 1. Original author :
559  !
560  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
561  !
562  ! 2. Last update :
563  !
564  ! See revisions.
565  !
566  ! 3. Revisions :
567  !
568  ! 20-Feb-2012 : Origination
569  !/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
570  !
571  ! 4. Copyright :
572  !
573  ! 5. Purpose :
574  !
575  ! Compute grid arrays for scrip for a specific unstructured grid
576  ! For interior vertices, we select for every triangle the barycenter
577  ! of the triangle. So to every node contained in N triangles we associate
578  ! a domain with N corners. Every one of those corners is contained
579  ! in 3 different domains.
580  ! For nodes that are on the boundary, we have to proceed differently.
581  ! For every such node, we have NEIGHBOR_PREV and NEIGHBOR_NEXT which
582  ! are the neighbor on each side of the boundary.
583  ! We put a corner on the middle of the edge. We also put a corner
584  ! on the node itself.
585  ! Note that instead of taking barycenters, we could have taken
586  ! Voronoi vertices, but this carries danger since Voronoi vertices
587  ! can be outside of the triangle. And it leaves open how to treat
588  ! the boundary.
589  !
590  ! 6. Method :
591  !
592  ! 7. Parameters, Variables and types :
593  !
594  ! 8. Called by :
595  !
596  ! Subroutine get_scrip_info
597  !
598  ! 9. Subroutines and functions used :
599  !
600  ! 10. Error messages:
601  !
602  ! 11. Remarks :
603  !
604  ! 12. Structure :
605  !
606  ! 13. Switches :
607  !
608  ! 14. Source code :
609  USE w3servmd, ONLY: extcde
610  USE w3gdatmd, ONLY : grids
611  IMPLICIT NONE
612  INTEGER, INTENT(IN) :: ID_GRD
613  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
614  real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
615  LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
616  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
617  real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
618  INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
619  INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
620 
621  INTEGER DIRAPPROACH, DUALAPPROACH, THEAPPROACH
622  INTEGER MNE, MNP, IE, IP, I
623  INTEGER NBPLUS, NBMINUS
624  INTEGER I1, I2, I3
625  real*8 :: elon1, elon2, elon3, elon, elonc
626  real*8 :: elat1, elat2, elat3, elat, elatc
627  REAL *8 :: DELTALON12, DELTALON13, DELTALAT12, DELTALAT13
628  REAL *8 :: THEDET
629  real*8 :: pt(3,2)
630  INTEGER, POINTER :: IOBP(:), TRIGINCD(:)
631  INTEGER, POINTER :: NEIGHBOR_PREV(:), NEIGHBOR_NEXT(:)
632  INTEGER, POINTER :: NBASSIGNEDCORNER(:), LISTNBCORNER(:)
633  INTEGER, POINTER :: STATUS(:), NEXTVERT(:), PREVVERT(:), FINALVERT(:)
634  INTEGER :: MAXCORNER, NBCORNER
635  INTEGER :: IDX, IPNEXT, IPPREV, NB, INEXT, IPREV
636  real*8, POINTER :: lon_cent_trig(:), lat_cent_trig(:)
637  real*8 :: elonip, elonnext, elonprev, elonn, elonp
638  real*8 :: elatip, elatnext, elatprev, elatn, elatp
639  INTEGER :: ISFINISHED, ZPREV
640  INTEGER :: DODEBUG
641  grid_rank=2
642  dirapproach=1
643  dualapproach=2
644  theapproach=dualapproach
645  mne=grids(id_grd)%NTRI
646  mnp=grids(id_grd)%NX
647  IF (theapproach .EQ. dirapproach) THEN
648  ALLOCATE(grid_center_lon(mne), stat=istat)
649  check_alloc_status( istat )
650  ALLOCATE(grid_center_lat(mne), stat=istat)
651  check_alloc_status( istat )
652  ALLOCATE(grid_corner_lon(3,mne), stat=istat)
653  check_alloc_status( istat )
654  ALLOCATE(grid_corner_lat(3,mne), stat=istat)
655  check_alloc_status( istat )
656  ALLOCATE(grid_mask(mne), stat=istat)
657  check_alloc_status( istat )
658  DO ie=1,mne
659  i1=grids(id_grd)%TRIGP(1,ie)
660  i2=grids(id_grd)%TRIGP(2,ie)
661  i3=grids(id_grd)%TRIGP(3,ie)
662  elon1=grids(id_grd)%XGRD(1,i1)
663  elon2=grids(id_grd)%XGRD(1,i2)
664  elon3=grids(id_grd)%XGRD(1,i3)
665  elat1=grids(id_grd)%YGRD(1,i1)
666  elat2=grids(id_grd)%YGRD(1,i2)
667  elat3=grids(id_grd)%YGRD(1,i3)
668  elon=(elon1 + elon2 + elon3)/3
669  elat=(elat1 + elat2 + elat3)/3
670  grid_center_lon(ie)=elon
671  grid_center_lat(ie)=elat
672  grid_corner_lon(1,ie)=elon1
673  grid_corner_lon(2,ie)=elon2
674  grid_corner_lon(3,ie)=elon3
675  grid_corner_lat(1,ie)=elat1
676  grid_corner_lat(2,ie)=elat2
677  grid_corner_lat(3,ie)=elat3
678  grid_mask(ie)=.true.
679  END DO
680  grid_corners=3
681  END IF
682  IF (theapproach .EQ. dualapproach) THEN
683  ALLOCATE(trigincd(mnp), stat=istat)
684  check_alloc_status( istat )
685  ALLOCATE(iobp(mnp), stat=istat)
686  check_alloc_status( istat )
687  ALLOCATE(neighbor_next(mnp), stat=istat)
688  check_alloc_status( istat )
689  ALLOCATE(neighbor_prev(mnp), stat=istat)
690  check_alloc_status( istat )
691  ALLOCATE(nbassignedcorner(mnp), stat=istat)
692  check_alloc_status( istat )
693  ALLOCATE(listnbcorner(mnp), stat=istat)
694  check_alloc_status( istat )
695 
696  ALLOCATE(status(mnp), stat=istat)
697  check_alloc_status( istat )
698  ALLOCATE(nextvert(mnp), stat=istat)
699  check_alloc_status( istat )
700  ALLOCATE(prevvert(mnp), stat=istat)
701  check_alloc_status( istat )
702  ALLOCATE(finalvert(mnp), stat=istat)
703  check_alloc_status( istat )
704  ALLOCATE(lon_cent_trig(mne), stat=istat)
705  check_alloc_status( istat )
706  ALLOCATE(lat_cent_trig(mne), stat=istat)
707  check_alloc_status( istat )
708 
709  CALL get_unstructured_vertex_degree (mnp, mne, &
710  grids(id_grd)%TRIGP, trigincd)
711  CALL get_boundary(mnp, mne, grids(id_grd)%TRIGP, iobp, &
712  neighbor_prev, neighbor_next)
713 
714  ! Find max number of corners
715  maxcorner=0
716  DO ip=1,mnp
717  IF (neighbor_next(ip) .EQ. 0) THEN
718  nbcorner=trigincd(ip)
719  ELSE
720  nbcorner=trigincd(ip) + 3
721  END IF
722  listnbcorner(ip)=nbcorner
723  IF (nbcorner .GT. maxcorner) THEN
724  maxcorner=nbcorner
725  END IF
726  END DO
727  grid_corners=maxcorner
728 
729  ALLOCATE(grid_center_lon(mnp), stat=istat)
730  check_alloc_status( istat )
731  ALLOCATE(grid_center_lat(mnp), stat=istat)
732  check_alloc_status( istat )
733  ALLOCATE(grid_corner_lon(maxcorner,mnp), stat=istat)
734  check_alloc_status( istat )
735  ALLOCATE(grid_corner_lat(maxcorner,mnp), stat=istat)
736  check_alloc_status( istat )
737  ALLOCATE(grid_mask(mnp), stat=istat)
738  check_alloc_status( istat )
739 
740  ! Add first three corners for boundaries
741  nbassignedcorner(:)=0
742  DO ip=1,mnp
743  grid_mask(ip)=.true.
744  IF (neighbor_next(ip) .GT. 0) THEN
745  ipnext=neighbor_next(ip)
746  ipprev=neighbor_prev(ip)
747  elonip=dble(grids(id_grd)%XGRD(1,ip))
748  elatip=dble(grids(id_grd)%YGRD(1,ip))
749  elonnext=dble(grids(id_grd)%XGRD(1,ipnext))
750  elatnext=dble(grids(id_grd)%YGRD(1,ipnext))
751  elonprev=dble(grids(id_grd)%XGRD(1,ipprev))
752  elatprev=dble(grids(id_grd)%YGRD(1,ipprev))
753 
754  ! Periodicity fix for corner node
755  IF ( abs(elonip - elonnext) .GT. 180.0 ) THEN
756  elonnext = elonnext -sign(360.0d0,(elonip - elonnext))
757  ENDIF
758  IF ( abs(elonip - elonprev) .GT. 180.0 ) THEN
759  elonprev = elonprev -sign(360.0d0,(elonip - elonprev))
760  ENDIF
761 
762  elonn=(elonip+elonnext)/2.0
763  elatn=(elatip+elatnext)/2.0
764  elonp=(elonip+elonprev)/2.0
765  elatp=(elatip+elatprev)/2.0
766 
767 
768  grid_corner_lon(1,ip)=elonn
769  grid_corner_lat(1,ip)=elatn
770  grid_corner_lon(2,ip)=elonip
771  grid_corner_lat(2,ip)=elatip
772  grid_corner_lon(3,ip)=elonp
773  grid_corner_lat(3,ip)=elatp
774  nbassignedcorner(ip)=3
775  END IF
776  END DO
777 
778  ! Compute centers
779  DO ip=1,mnp
780  grid_center_lon(ip)=dble(grids(id_grd)%XGRD(1,ip))
781  grid_center_lat(ip)=dble(grids(id_grd)%YGRD(1,ip))
782  END DO
783 
784  ! Check triangle node orientation
785  ! Compute triangle centers
786  nbplus=0
787  nbminus=0
788  DO ie=1,mne
789  i1=grids(id_grd)%TRIGP(1,ie)
790  i2=grids(id_grd)%TRIGP(2,ie)
791  i3=grids(id_grd)%TRIGP(3,ie)
792  pt(1,1)=dble(grids(id_grd)%XGRD(1,i1))
793  pt(2,1)=dble(grids(id_grd)%XGRD(1,i2))
794  pt(3,1)=dble(grids(id_grd)%XGRD(1,i3))
795  pt(1,2)=dble(grids(id_grd)%YGRD(1,i1))
796  pt(2,2)=dble(grids(id_grd)%YGRD(1,i2))
797  pt(3,2)=dble(grids(id_grd)%YGRD(1,i3))
798 
799  CALL fix_periodcity(pt)
800 
801  elon1 = pt(1,1)
802  elon2 = pt(2,1)
803  elon3 = pt(3,1)
804  elat1 = pt(1,2)
805  elat2 = pt(2,2)
806  elat3 = pt(3,2)
807 
808  deltalon12=elon2 - elon1
809  deltalon13=elon3 - elon1
810  deltalat12=elat2 - elat1
811  deltalat13=elat3 - elat1
812  thedet=deltalon12*deltalat13 - deltalon13*deltalat12
813  IF (thedet.GT.0) THEN
814  nbplus=nbplus+1
815  END IF
816  IF (thedet.LT.0) THEN
817  nbminus=nbminus+1
818  END IF
819  elon=(elon1 + elon2 + elon3)/3.0
820  elat=(elat1 + elat2 + elat3)/3.0
821 
822 
823  lon_cent_trig(ie)=elon
824  lat_cent_trig(ie)=elat
825 
826  END DO
827  dodebug=0
828  IF (dodebug.EQ.1) THEN
829  print *, 'nbplus=', nbplus, ' nbminus=', nbminus
830  END IF
831 
832  status(:) = 0
833  nextvert(:) = 0
834  prevvert(:) = 0
835  DO ie=1,mne
836  DO i=1,3
837  CALL triang_indexes(i, inext, iprev)
838  ip=grids(id_grd)%TRIGP(i,ie)
839  ipnext=grids(id_grd)%TRIGP(inext,ie)
840  ipprev=grids(id_grd)%TRIGP(iprev,ie)
841  IF (status(ip).EQ.0) THEN
842  IF (neighbor_next(ip).EQ.0) THEN
843  status(ip)=1
844  finalvert(ip)=ipprev
845  prevvert(ip)=ipprev
846  nextvert(ip)=ipnext
847  ELSE
848  IF (neighbor_prev(ip).EQ.ipprev) THEN
849  status(ip)=1
850  prevvert(ip)=ipprev
851  nextvert(ip)=ipnext
852  finalvert(ip)=neighbor_next(ip)
853  END IF
854  END IF
855  END IF
856  END DO
857  END DO
858  status(:)=0
859  DO
860  isfinished=1
861  DO ie=1,mne
862  elon=lon_cent_trig(ie)
863  elat=lat_cent_trig(ie)
864  DO i=1,3
865  CALL triang_indexes(i, inext, iprev)
866  ip=grids(id_grd)%TRIGP(i,ie)
867  ipnext=grids(id_grd)%TRIGP(inext,ie)
868  ipprev=grids(id_grd)%TRIGP(iprev,ie)
869  IF (status(ip).EQ.0) THEN
870  isfinished=0
871  zprev=prevvert(ip)
872  IF (zprev.EQ.ipprev) THEN
873  idx=nbassignedcorner(ip)
874  idx=idx+1
875  grid_corner_lon(idx,ip)=elon
876  grid_corner_lat(idx,ip)=elat
877  nbassignedcorner(ip)=idx
878  prevvert(ip)=ipnext
879  IF (ipnext.EQ.finalvert(ip)) THEN
880  status(ip)=1
881  END IF
882  END IF
883  END IF
884  END DO
885  END DO
886  IF (isfinished.EQ.1) THEN
887  EXIT
888  END IF
889  END DO
890  DO ip=1,mnp
891  IF (nbassignedcorner(ip).NE.listnbcorner(ip)) THEN
892  WRITE(*,*) 'Incoherent number at IP=', ip
893  WRITE(*,*) ' NbAssignedCorner(IP)=', nbassignedcorner(ip)
894  WRITE(*,*) ' ListNbCorner(IP)=', listnbcorner(ip)
895  WRITE(*,*) ' N_N=', neighbor_next(ip), 'N_P=', neighbor_prev(ip)
896  WRITE(*,*) ' TrigIncd=', trigincd(ip)
897  stop 'wmscrpmd, case 2'
898  END IF
899  END DO
900 
901  ! if the number of corner is below threshold, we have to
902  ! add some more.
903  DO ip=1,mnp
904  nb=nbassignedcorner(ip)
905  IF (nb .LT. maxcorner) THEN
906  elon=grid_corner_lon(nb,ip)
907  elat=grid_corner_lat(nb,ip)
908  DO idx=nb+1,maxcorner
909  grid_corner_lon(idx,ip)=elon
910  grid_corner_lat(idx,ip)=elat
911  END DO
912  END IF
913  END DO
914  DEALLOCATE(nbassignedcorner, stat=istat)
915  check_dealloc_status( istat )
916  DEALLOCATE(listnbcorner, stat=istat)
917  check_dealloc_status( istat )
918  DEALLOCATE(trigincd, stat=istat)
919  check_dealloc_status( istat )
920  DEALLOCATE(iobp, stat=istat)
921  check_dealloc_status( istat )
922  DEALLOCATE(neighbor_prev, stat=istat)
923  check_dealloc_status( istat )
924  DEALLOCATE(neighbor_next, stat=istat)
925  check_dealloc_status( istat )
926  DEALLOCATE(status, stat=istat)
927  check_dealloc_status( istat )
928  DEALLOCATE(nextvert, stat=istat)
929  check_dealloc_status( istat )
930  DEALLOCATE(prevvert, stat=istat)
931  check_dealloc_status( istat )
932  DEALLOCATE(finalvert, stat=istat)
933  check_dealloc_status( istat )
934  DEALLOCATE(lon_cent_trig, stat=istat)
935  check_dealloc_status( istat )
936  DEALLOCATE(lat_cent_trig, stat=istat)
937  check_dealloc_status( istat )
938 
939  ALLOCATE(grid_dims(2), stat=istat)
940  check_alloc_status( istat )
941  grid_dims(1)=mnp
942  grid_dims(2)=1
943  grid_size=mnp
944  END IF

References w3servmd::extcde(), fix_periodcity(), get_boundary(), get_unstructured_vertex_degree(), w3gdatmd::grids, and triang_indexes().

Referenced by get_scrip_info().

◆ get_unstructured_vertex_degree()

subroutine wmscrpmd::get_unstructured_vertex_degree ( integer, intent(in)  MNP,
integer, intent(in)  MNE,
integer, dimension(:,:), intent(in)  TRIGP,
integer, dimension(:), intent(out)  TRIGINCD 
)

This function returns the list of incidences.

Output: TrigIncd - number of triangles contained by vertices.

Parameters
[in]MNPNumber of nodes
[in]MNEList of nodes
[in]TRIGPNumber of triangles
[out]TRIGINCDNumber of triangles contained by vertices.
Author
M. Dutour
A. Roland
Date
20-Feb-2012

Definition at line 1448 of file wmscrpmd.F90.

1448  ! Written:
1449  !
1450  ! 20-Feb-2012
1451  !
1452  ! Author:
1453  !
1454  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
1455  !
1456  ! Parameters:
1457  ! Input:
1458  ! MNP: number of nodes
1459  ! INE: list of nodes
1460  ! MNE: number of triangles
1461  ! Output:
1462  ! TrigIncd (number of triangles contained by vertices
1463  !
1464  ! Description:
1465  ! this function returns the list of incidences
1466  !
1467  IMPLICIT NONE
1468  INTEGER, INTENT(IN) :: MNP, MNE
1469  INTEGER, INTENT(IN) :: TRIGP(:,:)
1470  INTEGER, INTENT(OUT) :: TRIGINCD(:)
1471  INTEGER :: IP, IE, I
1472  trigincd=0
1473  DO ie=1,mne
1474  DO i=1,3
1475  ip=trigp(i,ie)
1476  trigincd(ip)=trigincd(ip) + 1
1477  END DO
1478  END DO

Referenced by get_scrip_info_unstructured().

◆ scrip_info_renormalization()

subroutine wmscrpmd::scrip_info_renormalization ( real*8, dimension(:), intent(inout)  GRID_CENTER_LON,
real*8, dimension(:), intent(inout)  GRID_CENTER_LAT,
real*8, dimension(:,:), intent(inout)  GRID_CORNER_LON,
real*8, dimension(:,:), intent(inout)  GRID_CORNER_LAT,
logical, dimension(:), intent(in)  GRID_MASK,
integer, dimension(:), intent(in)  GRID_DIMS,
integer, intent(in)  GRID_SIZE,
integer, intent(in)  GRID_CORNERS,
integer, intent(in)  GRID_RANK,
real*8  CONV_DX,
real*8  CONV_DY,
real*8  OFFSET,
real*8  GRIDSHIFT 
)

Rescale according to whether the grid is spherical or not.

This is adapted from Erick Rogers scrip_wrapper.

Purpose is to rescale according to whether the grid is spherical or not and to adjust by some small shift to keep SCRIP happy in situations where nodes of different grids overlay.

We apply various transformations to the longitude latitude following here the transformations that were done only in finite difference meshes.

Parameters
[in,out]GRID_CENTER_LON
[in,out]GRID_CENTER_LAT
[in,out]GRID_CORNER_LON
[in,out]GRID_CORNER_LAT
[in]GRID_MASK
[in]GRID_DIMS
[in]GRID_SIZE
[in]GRID_CORNERS
[in]GRID_RANK
CONV_DX
CONV_DY
OFFSET
GRIDSHIFT
Author
M. Dutour
A. Roland
Date
20-Feb-2012

Definition at line 1313 of file wmscrpmd.F90.

1313  ! 1. Original author :
1314  !
1315  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
1316  ! Adapted from Erick Rogers scrip_wrapper
1317  !
1318  ! 2. Last update :
1319  !
1320  ! See revisions.
1321  !
1322  ! 3. Revisions :
1323  !
1324  ! 20-Feb-2012 : Origination
1325  !
1326  ! 4. Copyright :
1327  !
1328  ! 5. Purpose :
1329  !
1330  ! This is adapted from Erick Rogers scrip_wrapper
1331  ! Purpose is to rescale according to whether the grid is spherical
1332  ! or not and to adjust by some small shift to keep SCRIP happy
1333  ! in situations where nodes of different grids overlay
1334  !
1335  ! 6. Method :
1336  !
1337  ! We apply various transformations to the longitude latitude
1338  ! following here the transformations that were done only in
1339  ! finite difference meshes.
1340  !
1341  ! 7. Parameters, Variables and types :
1342  !
1343  ! 8. Called by :
1344  !
1345  ! Subroutine WMGHGH
1346  !
1347  ! 9. Subroutines and functions used :
1348  !
1349  ! 10. Error messages:
1350  !
1351  ! 11. Remarks :
1352  !
1353  ! 12. Structure :
1354  !
1355  ! 13. Switches :
1356  !
1357  ! 14. Source code :
1358  IMPLICIT NONE
1359  real*8, INTENT(INOUT) :: grid_center_lon(:)
1360  real*8, INTENT(INOUT) :: grid_center_lat(:)
1361  LOGICAL, INTENT(IN) :: GRID_MASK(:)
1362  real*8, INTENT(INOUT) :: grid_corner_lon(:,:)
1363  real*8, INTENT(INOUT) :: grid_corner_lat(:,:)
1364  INTEGER, INTENT(IN) :: GRID_DIMS(:)
1365  INTEGER, INTENT(IN) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
1366  real*8 :: conv_dx, conv_dy, offset, gridshift
1367  real*8 deg2rad
1368  !
1369  INTEGER :: I, J, IP
1370  real*8 :: minlon, minlat, maxlon, maxlat, hlon, hlat
1371  real*8 :: minloncorner, maxloncorner, minlatcorner, maxlatcorner
1372 
1373  DO i=1,grid_size
1374  grid_center_lon(i)=(grid_center_lon(i)+offset)/conv_dx + &
1375  & gridshift
1376  grid_center_lat(i)=grid_center_lat(i)/conv_dy + &
1377  & gridshift
1378  IF(grid_center_lon(i)>360.0) THEN
1379  grid_center_lon(i)=grid_center_lon(i)-360.0
1380  END IF
1381  IF(grid_center_lon(i)<000.0) THEN
1382  grid_center_lon(i)=grid_center_lon(i)+360.0
1383  END IF
1384  DO j=1,grid_corners
1385  grid_corner_lon(j, i)=(grid_corner_lon(j, i)+offset)/conv_dx+ &
1386  & gridshift
1387  grid_corner_lat(j, i)=grid_corner_lat(j, i)/conv_dy + &
1388  & gridshift
1389  IF(grid_corner_lon(j,i)>360.0) THEN
1390  grid_corner_lon(j,i)=grid_corner_lon(j,i)-360.0
1391  END IF
1392  IF(grid_corner_lon(j,i)<000.0) THEN
1393  grid_corner_lon(j,i)=grid_corner_lon(j,i)+360.0
1394  END IF
1395  END DO
1396  END DO
1397 

Referenced by scrip_wrapper().

◆ scrip_wrapper()

subroutine wmscrpmd::scrip_wrapper ( integer(scrip_i4), intent(in)  ID_SRC,
integer(scrip_i4), intent(in)  ID_DST,
integer(scrip_i4), dimension(:,:), intent(in)  MAPSTA_SRC,
integer(scrip_i4), dimension(:,:), intent(in)  MAPST2_SRC,
logical(scrip_logical), intent(in)  FLAGLL,
real (scrip_r8), intent(in)  GRIDSHIFT,
logical(scrip_logical), intent(in)  L_MASTER,
logical(scrip_logical), intent(in)  L_READ,
logical(scrip_logical), intent(in)  L_TEST 
)

Compute grid information required by SCRIP.

Parameters
[in]ID_SRC
[in]ID_DST
[in]MAPSTA_SRC
[in]MAPST2_SRC
[in]FLAGLL
[in]GRIDSHIFT
[in]L_MASTER
[in]L_READ
[in]L_TEST
Author
E. Rogers
M. Dutour
A. Roland
Date
10-Dec-2014

Definition at line 115 of file wmscrpmd.F90.

115  !/ +-----------------------------------+
116  !/ | WAVEWATCH III |
117  !/ | E. Rogers, M. Dutour, A. Roland |
118  !/ | FORTRAN 90 |
119  !/ | Last update : 10-Dec-2014 !
120  !/ +-----------------------------------+
121  !/
122  ! 1. Original author :
123  !
124  ! Erick Rogers, NRL
125  !
126  ! 2. Last update :
127  !
128  ! See revisions.
129  !
130  ! 3. Revisions :
131  !
132  ! 29-Apr-2011 : Origination ( version 4.01 )
133  ! 20-Feb-2012 : Mathieu Dutour Sikiric, Aron Roland Z&P
134  ! Modification is to split the code into several
135  ! subroutines
136  ! ---get_scrip_info_structured (the structured case)
137  ! ---get_scrip_info (chooses according to FD/FE)
138  ! ---scrip_info_renormalization (conv_x and all that)
139  ! 11-Apr-2013 Kevin Lind
140  ! Added code for reading/writing SCRIP remap files
141  !/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
142  !
143  ! 4. Copyright :
144  !
145  ! 5. Purpose :
146  !
147  ! Compute grid information required by SCRIP
148  !
149  ! 6. Method :
150  !
151  ! 7. Parameters, Variables and types :
152  !
153  ! 8. Called by :
154  !
155  ! Subroutine WMGHGH
156  !
157  ! 9. Subroutines and functions used :
158  !
159  ! Subroutine SCRIP
160  !
161  ! 10. Error messages:
162  !
163  ! 11. Remarks :
164  !
165  ! 12. Structure :
166  !
167  ! 13. Switches :
168  !
169  ! 14. Source code :
170 
171  USE scrip_grids
172  USE scrip_remap_vars
173  USE scrip_constants
174  USE scrip_kindsmod
175  USE scrip_interface
176  USE w3servmd, ONLY: extcde
177  ! USE W3GDATMD, ONLY : GRIDS
178 
179  IMPLICIT NONE
180  INTEGER(SCRIP_I4), INTENT(IN) :: ID_SRC, ID_DST
181  INTEGER(SCRIP_I4), INTENT(IN) :: MAPSTA_SRC(:,:)
182  INTEGER(SCRIP_I4), INTENT(IN) :: MAPST2_SRC(:,:)
183  LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: FLAGLL
184  REAL (SCRIP_R8), INTENT(IN) :: GRIDSHIFT
185  LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_MASTER ! Am I the master processor (do I/O)?
186  LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_READ ! Do I read the remap file?
187  LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_TEST ! Whether to include test output
188  ! in subroutines
189  !/ ------------------------------------------------------------------- /
190  !/ local variables
191  !/
192  INTEGER(SCRIP_I4) :: IREC,I,J,NI,NJ,IDUM,NK,K, &
193  ILINK,IW,ICORNER,NGOODPNTS, &
194  NBADPNTS
195  INTEGER(SCRIP_I4) :: ISRC,JSRC,KSRC,IPNT,KDST, &
196  NI_SRC
197  REAL (SCRIP_R8) :: LAT_CONVERSION,OFFSET
198  REAL (SCRIP_R8) :: CONV_DX,CONV_DY,WEIGHT
199  REAL (SCRIP_R8) :: WTSUM
200 #ifdef W3_T38
201  CHARACTER (LEN=10) :: CDATE_TIME(3)
202  INTEGER :: DATE_TIME(8)
203  INTEGER :: ELAPSED_TIME, BEG_TIME, &
204  END_TIME
205 #endif
206 
207  ! test output for input variables
208 
209 #ifdef W3_T38
210  if(l_master)write(*,*)'flagll = ',flagll
211  if(l_master)write(*,*)'gridshift = ',gridshift
212 #endif
213 
214  !
215  ! START: universal settings
216  !
217 
218  ! Set variables for converting to degrees
219  ! notes: SCRIP only operates on spherical coordinates, so for the case
220  ! where the problem is specified by the user as in a
221  ! meters/cartesian coordinate system, it is necessary to make
222  ! a temporary conversion to a "fake" spherical coordinate grid,
223  ! to keep SCRIP happy. The good news here is that multi-grid
224  ! meters-grid simulations will be very rare: we will probably only
225  ! encounter them in the context of simple test cases. Strictly
226  ! speaking, this conversion does not even need to be physically
227  ! correct, e.g. we could say that 1000 km is 1 deg....as long as
228  ! we are consistent between grids.
229  ! Potential future improvement: make conv_dy and offset such that dst grid
230  ! covers a specific longitude range, e.g. 1 deg east to 2 deg east
231 
232  !
233  ! START: set up src grid
234  !
235 
236  !notes: when we work out how to interface with an unstructured grid,
237  ! we will need to revisit this issue of how to set grid1_rank, etc.
238  ! strategy: declare variables in grid module, but allocate them here.
239  !
240  grid1_units='degrees' ! the other option is radians...we don't use this
241  grid1_name='src' ! this is not used, except for netcdf output
242  CALL get_scrip_info(id_src, &
243  & grid1_center_lon, grid1_center_lat, &
244  & grid1_corner_lon, grid1_corner_lat, grid1_mask, &
245  & grid1_dims, grid1_size, grid1_corners, grid1_rank)
246  grid2_units='degrees'
247  grid2_name='dst'
248  CALL get_scrip_info(id_dst, &
249  & grid2_center_lon, grid2_center_lat, &
250  & grid2_corner_lon, grid2_corner_lat, grid2_mask, &
251  & grid2_dims, grid2_size, grid2_corners, grid2_rank)
252 
253  IF(flagll)THEN
254  conv_dx=one
255  conv_dy=one
256  offset=zero
257  ELSE
258  lat_conversion=zero ! lat_conversion
259  ! notes: this is the latitude used for conversion everywhere
260  ! in the grid (approximation)
261  ! (in radians)
262  ! conv_dy=92.6*1200.0 ! physical, =92.6/(3/3600)=111000 m = 111 km
263  conv_dy=1.0e+6_scrip_r8 ! non-physical, 1e+6=1 deg
264  conv_dx=cos(lat_conversion)*conv_dy
265  ! notes: offset (in meters), is necessary so that our grid does
266  ! not lie on the branch cut
267  offset=75000.0_scrip_r8-min(minval(grid1_center_lon), &
268  & minval(grid2_center_lon))
269  ENDIF
270 
271  !.....test output
272 #ifdef W3_T38
273  write(*,*)'l_master = ',l_master
274  if(l_master)then
275  write(*,*)'conv_dx=', conv_dx
276  write(*,*)'conv_dy=', conv_dy
277  write(*,*)'offset = ',offset
278  write(*,*)'grid1_size=', grid1_size
279  write(*,*)'grid2_size=', grid2_size
280  write(*,*)'l_read = ',l_read
281  write(*,*)'minval(grid1_center_lon) = ',minval(grid1_center_lon)
282  write(*,*)'maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
283  write(*,*)'minval(grid1_center_lat) = ',minval(grid1_center_lat)
284  write(*,*)'maxval(grid1_center_lat) = ',maxval(grid1_center_lat)
285  write(*,*)'minval(grid2_center_lon) = ',minval(grid2_center_lon)
286  write(*,*)'maxval(grid2_center_lon) = ',maxval(grid2_center_lon)
287  write(*,*)'minval(grid2_center_lat) = ',minval(grid2_center_lat)
288  write(*,*)'maxval(grid2_center_lat) = ',maxval(grid2_center_lat)
289  endif
290 #endif
291 
292  CALL scrip_info_renormalization( &
293  & grid1_center_lon, grid1_center_lat, &
294  & grid1_corner_lon, grid1_corner_lat, grid1_mask, &
295  & grid1_dims, grid1_size, grid1_corners, grid1_rank, &
296  & conv_dx, conv_dy, offset, gridshift)
297  CALL scrip_info_renormalization( &
298  & grid2_center_lon, grid2_center_lat, &
299  & grid2_corner_lon, grid2_corner_lat, grid2_mask, &
300  & grid2_dims, grid2_size, grid2_corners, grid2_rank, &
301  & conv_dx, conv_dy, offset, zero)
302 
303  !.....Set constants for thresholding weights:
304  frac_lowest =1.e-3_scrip_r8
305  frac_highest=one+1.e-3_scrip_r8
306  wt_lowest =zero
307  wt_highest =one+1.e-7_scrip_r8
308 
309 #ifdef W3_T38
310  call date_and_time (cdate_time(1), cdate_time(2), cdate_time(3), date_time)
311  beg_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
312 #endif
313  CALL scrip(id_src, id_dst, l_master, l_read, l_test)
314 #ifdef W3_T38
315  call date_and_time (cdate_time(1), cdate_time(2), cdate_time(3), date_time)
316  end_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
317  elapsed_time = end_time - beg_time
318  write(0,*) "SCRIP: ", elapsed_time, " MSEC"
319 #endif
320 
321 #ifdef W3_T38
322  if(l_master)write(*,*)'new minval(grid1_center_lon) = ',minval(grid1_center_lon)
323  if(l_master)write(*,*)'new maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
324 #endif
325 
326  !.....notes: at this point we have the following useful variables:
327  ! num_wts, e.g. num_wts=3....for first order conservative remapping,
328  ! only the first one is relevant, the other two are for second-order
329  ! remapping.
330  ! max_links_map1, e.g. max_links_map1=1369,
331  ! grid2_size, e.g. grid2_size=1849,
332  ! wts_map1(num_wts,max_links_map1), e.g. wts_map1(3,1369),
333  ! grid1_add_map1(max_links_map1), e.g. grid1_add_map1(1369),
334  ! grid2_add_map1(max_links_map1), e.g. grid2_add_map1(1369),
335  ! grid2_frac(grid2_size), e.g. grid2_frac(1849),
336  ! (see earlier versions for notes re: equivalency in netcdf/matlab)
337  !
338  !.....test output (optional)
339  !
340  !.....note re: notation: I use k for the combined i/j array, similar to isea,
341  ! but not necessarily the same as isea since some points may
342  ! be land etc.
343 #ifdef W3_T38
344  if(l_master)then
345  do k=1,grid2_size
346  write(403,*)grid2_frac(k)
347  end do
348  do ilink=1,max_links_map1
349  write(405,'(999(1x,f20.7))')(wts_map1(iw,ilink),iw=1,num_wts)
350  end do
351  do ilink=1,max_links_map1
352  write(406,'(i20)')grid1_add_map1(ilink) ! equivalent to
353  ! my "src_address"
354  write(407,'(i20)')grid2_add_map1(ilink) ! equivalent to
355  ! my "dst_address"
356  end do
357  endif
358 #endif
359 
360  !.....organize results and return to wmghgh.
361 
362  ! what we need, for purpose of feeding back to ww3, for each dst grid node:
363  ! a) the set of src grid nodes, in terms of isea, for which
364  ! weights are available
365  ! b) the corresponding set of weights
366 
367  ALLOCATE ( wgtdata(grid2_size), stat=istat )
368  check_alloc_status( istat )
369 
370  !.....step 1: count up NR0, NR1, NR2, NRL, NLOC (NR1 and NLOC are denoted
371  ! "n" here)
372  ! It is especially important to determine how large npnts gets,
373  ! so that we can allocate arrays properly
374  wgtdata%NR0=0
375  wgtdata%NR2=0
376  wgtdata%NRL=0
377  wgtdata%N=0
378 
379  ni_src=grid1_dims(1)
380  ngoodpnts=0
381  nbadpnts=0
382 
383  DO ilink=1,max_links_map1
384 
385  !........note that this pair of if-thens *must* be consistent with the
386  !........single if-then below
387  IF ((grid2_frac(grid2_add_map1(ilink))>frac_lowest) .AND. &
388  (grid2_frac(grid2_add_map1(ilink))<frac_highest)) THEN
389  ! then consider this link either to include, or to print a warning
390  IF( (wts_map1(1,ilink)>=wt_lowest) .AND. &
391  (wts_map1(1,ilink)<=wt_highest) ) THEN
392  ksrc=grid1_add_map1(ilink)
393  jsrc=int((ksrc-1)/ni_src)+1
394  isrc=ksrc-(jsrc-1)*ni_src
395 
396  IF (mapsta_src(jsrc,isrc).EQ.0) THEN ! excluded point
397  wgtdata(grid2_add_map1(ilink))%NR0 = &
398  wgtdata(grid2_add_map1(ilink))%NR0 + 1
399  IF (mapst2_src(jsrc,isrc).EQ.0)THEN
400  wgtdata(grid2_add_map1(ilink))%NRL = &
401  wgtdata(grid2_add_map1(ilink))%NRL + 1
402  ENDIF
403  ELSE IF (abs(mapsta_src(jsrc,isrc)).EQ.1) THEN
404  ! sea point
405  wgtdata(grid2_add_map1(ilink))%N= &
406  wgtdata(grid2_add_map1(ilink))%N+1
407  ELSE IF (abs(mapsta_src(jsrc,isrc)).EQ.2) THEN
408  ! bnd point
409  wgtdata(grid2_add_map1(ilink))%NR2 = &
410  wgtdata(grid2_add_map1(ilink))%NR2 + 1
411  END IF
412  ngoodpnts=ngoodpnts+1
413  ELSEIF ( (grid1_frac(grid1_add_map1(ilink))>frac_lowest) &
414  .AND. (grid1_frac(grid1_add_map1(ilink))<frac_highest)&
415  ) THEN
416  ! note that we don't bother giving a warning for the
417  ! cases where grid1_frac is strange.
418  nbadpnts=nbadpnts+1
419  ! we also don't bother giving a warning if we've already
420  ! given a lot of warnings (we keep counting, though)
421  IF((nbadpnts.LT.5).AND.l_master)THEN ! print warnings
422  WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
423  WRITE(*,'(A)')'WARNING: SCRIP weight problem '
424  WRITE(*,'(4x,A,I7,A,I7)')'ilink = ',ilink,' out of ',&
425  max_links_map1
426  WRITE(*,'(4x,A,I7)')'grid1_add_map1(ilink) = ', &
427  grid1_add_map1(ilink)
428  WRITE(*,'(4x,A,I7)')'grid2_add_map1(ilink) = ', &
429  grid2_add_map1(ilink)
430  WRITE(*,'(4x,A,E12.4)')'wts_map1(1,ilink) = ', &
431  wts_map1(1,ilink)
432  WRITE(*,'(4x,A,F10.5)') &
433  'grid1_frac(grid1_add_map1(ilink)) = ', &
434  grid1_frac(grid1_add_map1(ilink))
435  WRITE(*,'(4x,A,F10.5)') &
436  'grid2_frac(grid2_add_map1(ilink)) = ', &
437  grid2_frac(grid2_add_map1(ilink))
438  WRITE(*,'(4x,A,F10.5)')'grid1_center_lat = ', &
439  grid1_center_lat(grid1_add_map1(ilink))
440  WRITE(*,'(4x,A,F10.5)')'grid1_center_lon = ', &
441  grid1_center_lon(grid1_add_map1(ilink))
442  WRITE(*,'(4x,A,F10.5)')'grid2_center_lat = ', &
443  grid2_center_lat(grid2_add_map1(ilink))
444  WRITE(*,'(4x,A,F10.5)')'grid2_center_lon = ', &
445  grid2_center_lon(grid2_add_map1(ilink))
446  WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
447  ENDIF
448  ENDIF
449  ENDIF
450  END DO
451  IF((nbadpnts.GT.0).AND.l_master)THEN
452  WRITE(*,'(4x,A,I5,A)')'We had problems in ',nbadpnts, &
453  ' points.'
454  WRITE(*,'(4x,I8,A)')ngoodpnts,' points appear to be OK.'
455  ENDIF
456  IF( (nbadpnts.GT.(ngoodpnts/30)) .AND.l_master )THEN
457  WRITE(*,'(4x,A)')'Error: excessive SCRIP failure. Stopping.'
458  stop 'wmscrpmd, case 1'
459  ENDIF
460 
461  !.....step 2: allocate according to the size "n" determined above
462  DO kdst=1,grid2_size
463  ALLOCATE ( wgtdata(kdst)%W(wgtdata(kdst)%N), stat=istat )
464  check_alloc_status( istat )
465  ALLOCATE ( wgtdata(kdst)%K(wgtdata(kdst)%N), stat=istat )
466  check_alloc_status( istat )
467  wgtdata(kdst)%N=0
468  END DO
469 
470  !.....step 3: save weights
471  DO ilink=1,max_links_map1
472 
473  ksrc=grid1_add_map1(ilink)
474  jsrc=int((ksrc-1)/ni_src)+1
475  isrc=ksrc-(jsrc-1)*ni_src
476 
477  !........note that this single if-then *must* be consistent with the
478  !........pair of if-thens above
479  IF ((grid2_frac(grid2_add_map1(ilink))>frac_lowest) .AND. &
480  & (grid2_frac(grid2_add_map1(ilink))<frac_highest) .AND. &
481  & (wts_map1(1,ilink)>=wt_lowest) .AND. &
482  & (wts_map1(1,ilink)<=wt_highest))THEN
483  IF (abs(mapsta_src(jsrc,isrc)).EQ.1) THEN ! sea point
484  wgtdata(grid2_add_map1(ilink))%N= &
485  & wgtdata(grid2_add_map1(ilink))%N+1
486  wgtdata(grid2_add_map1(ilink))%W(wgtdata( &
487  & grid2_add_map1(ilink))%N)=wts_map1(1,ilink)
488  wgtdata(grid2_add_map1(ilink))%K(wgtdata( &
489  & grid2_add_map1(ilink))%N)=grid1_add_map1(ilink)
490  ENDIF
491  ENDIF
492  END DO
493 
494  !.....step 4: re-normalize weights. This is necessary because we called
495  !.....scrip without the mask. Now that we have the mask in place, we need
496  !.....to re-normalize the weights.
497  DO kdst=1,grid2_size
498  IF (wgtdata(kdst)%N > 0) THEN
499  wtsum=zero
500  DO ipnt=1,wgtdata(kdst)%N
501  wtsum=wtsum+wgtdata(kdst)%W(ipnt)
502  END DO
503  DO ipnt=1,wgtdata(kdst)%N
504  wgtdata(kdst)%W(ipnt)=wgtdata(kdst)%W(ipnt)/wtsum
505  END DO
506  END IF
507  END DO
508 
509  CALL scrip_clear

References w3servmd::extcde(), scrip_remap_vars::frac_highest, scrip_remap_vars::frac_lowest, get_scrip_info(), scrip_remap_vars::grid1_add_map1, scrip_grids::grid1_center_lat, scrip_grids::grid1_center_lon, scrip_grids::grid1_corner_lat, scrip_grids::grid1_corner_lon, scrip_grids::grid1_corners, scrip_grids::grid1_dims, scrip_grids::grid1_frac, scrip_grids::grid1_mask, scrip_grids::grid1_name, scrip_grids::grid1_rank, scrip_grids::grid1_size, scrip_grids::grid1_units, scrip_remap_vars::grid2_add_map1, scrip_grids::grid2_center_lat, scrip_grids::grid2_center_lon, scrip_grids::grid2_corner_lat, scrip_grids::grid2_corner_lon, scrip_grids::grid2_corners, scrip_grids::grid2_dims, scrip_grids::grid2_frac, scrip_grids::grid2_mask, scrip_grids::grid2_name, scrip_grids::grid2_rank, scrip_grids::grid2_size, scrip_grids::grid2_units, scrip_remap_vars::max_links_map1, scrip_remap_vars::num_wts, scrip_constants::one, scrip_interface::scrip(), scrip_interface::scrip_clear(), scrip_info_renormalization(), scrip_interface::wgtdata, scrip_remap_vars::wt_highest, scrip_remap_vars::wt_lowest, scrip_remap_vars::wts_map1, and scrip_constants::zero.

Referenced by wmgridmd::wmghgh().

◆ triang_indexes()

subroutine wmscrpmd::triang_indexes ( integer, intent(in)  I,
integer, intent(out)  INEXT,
integer, intent(out)  IPREV 
)

Desc not available.

Parameters
[in]I
[out]INEXT
[out]IPREV
Author
M. Dutour
A. Roland
Date
NA

Definition at line 1413 of file wmscrpmd.F90.

1413  ! 1. Original author :
1414  !
1415  ! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
1416  !
1417  INTEGER, INTENT(IN) :: I
1418  INTEGER, INTENT(OUT) :: INEXT, IPREV
1419  IF (i.EQ.1) THEN
1420  inext=3
1421  ELSE
1422  inext=i-1
1423  END IF
1424  IF (i.EQ.3) THEN
1425  iprev=1
1426  ELSE
1427  iprev=i+1
1428  END IF

Referenced by get_boundary(), and get_scrip_info_unstructured().

scrip_grids::grid1_corners
integer(scrip_i4), save grid1_corners
Definition: scrip_grids.f:68
scrip_remap_vars::grid1_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid1_add_map1
Definition: scrip_remap_vars.f:77
w3triamd::triang_indexes
subroutine triang_indexes(I, INEXT, IPREV)
Set indices of the triangle.
Definition: w3triamd.F90:2720
scrip_constants::half
real(kind=scrip_r8), parameter half
Definition: scrip_constants.f:50
scrip_grids::grid1_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid1_mask
Definition: scrip_grids.f:93
scrip_remap_vars::wt_lowest
real(kind=scrip_r8) wt_lowest
Definition: scrip_remap_vars.f:88
w3gdatmd::ungtype
integer, parameter ungtype
Definition: w3gdatmd.F90:626
scrip_remap_vars::wt_highest
real(kind=scrip_r8) wt_highest
Definition: scrip_remap_vars.f:88
scrip_grids::grid2_units
character(scrip_charlength), save grid2_units
Definition: scrip_grids.f:80
scrip_grids::grid2_rank
integer(scrip_i4), save grid2_rank
Definition: scrip_grids.f:68
scrip_grids::grid2_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid2_mask
Definition: scrip_grids.f:93
scrip_grids::grid2_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lat
Definition: scrip_grids.f:103
w3gdatmd::grids
type(grid), dimension(:), allocatable, target grids
Definition: w3gdatmd.F90:1088
scrip_grids::deg2rad
real(scrip_r8), parameter deg2rad
Definition: scrip_grids.f:84
scrip_grids::grid1_rank
integer(scrip_i4), save grid1_rank
Definition: scrip_grids.f:68
scrip_constants::one
real(kind=scrip_r8), parameter one
Definition: scrip_constants.f:50
scrip_grids
Definition: scrip_grids.f:49
scrip_timers::status
character(len=8), dimension(max_timers), save status
Definition: scrip_timers.f:63
scrip_grids::grid2_corners
integer(scrip_i4), save grid2_corners
Definition: scrip_grids.f:68
wmscrpmd::scrip_info_renormalization
subroutine scrip_info_renormalization(GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK, CONV_DX, CONV_DY, OFFSET, GRIDSHIFT)
Rescale according to whether the grid is spherical or not.
Definition: wmscrpmd.F90:1313
scrip_constants::zero
real(kind=scrip_r8), parameter zero
Definition: scrip_constants.f:50
w3triamd::fix_periodcity
subroutine fix_periodcity(I1, I2, I3, XGRD, YGRD, PT)
Adjust element longitude coordinates for elements straddling the dateline with distance of ~360 degre...
Definition: w3triamd.F90:3059
scrip_interface
Definition: scrip_interface.F90:2
w3servmd
Definition: w3servmd.F90:3
scrip_grids::grid2_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lat
Definition: scrip_grids.f:120
scrip_grids::grid1_dims
integer(scrip_i4), dimension(:), allocatable, save grid1_dims
Definition: scrip_grids.f:74
scrip_remap_vars::max_links_map1
integer(scrip_i4), save max_links_map1
Definition: scrip_remap_vars.f:66
scrip_grids::grid1_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lon
Definition: scrip_grids.f:103
scrip_grids::grid1_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lat
Definition: scrip_grids.f:120
scrip_remap_vars::frac_highest
real(kind=scrip_r8) frac_highest
Definition: scrip_remap_vars.f:87
scrip_grids::grid2_size
integer(scrip_i4), save grid2_size
Definition: scrip_grids.f:68
scrip_grids::grid1_name
character(scrip_charlength), save grid1_name
Definition: scrip_grids.f:77
scrip_remap_vars::frac_lowest
real(kind=scrip_r8) frac_lowest
Definition: scrip_remap_vars.f:87
scrip_interface::scrip
subroutine scrip(src_num, dst_num, l_master, l_read, l_test)
Definition: scrip_interface.F90:330
scrip_grids::grid2_name
character(scrip_charlength), save grid2_name
Definition: scrip_grids.f:77
scrip_interface::wgtdata
type(weight_data), dimension(:), allocatable wgtdata
Definition: scrip_interface.F90:73
scrip_remap_vars::wts_map1
real(scrip_r8), dimension(:,:), allocatable, save wts_map1
Definition: scrip_remap_vars.f:83
scrip_constants
Definition: scrip_constants.f:38
scrip_grids::grid2_dims
integer(scrip_i4), dimension(:), allocatable, save grid2_dims
Definition: scrip_grids.f:74
scrip_grids::grid1_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lat
Definition: scrip_grids.f:103
scrip_kindsmod
Definition: scrip_kindsmod.f90:3
scrip_grids::grid1_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lon
Definition: scrip_grids.f:120
scrip_remap_vars::grid2_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid2_add_map1
Definition: scrip_remap_vars.f:77
wmscrpmd::get_scrip_info_structured
subroutine get_scrip_info_structured(ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
Compute grid arrays needed for scrip for a specific structured grid.
Definition: wmscrpmd.F90:974
w3gdatmd
Definition: w3gdatmd.F90:16
scrip_grids::grid1_size
integer(scrip_i4), save grid1_size
Definition: scrip_grids.f:68
wmscrpmd::get_scrip_info_unstructured
subroutine get_scrip_info_unstructured(ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
Compute grid arrays for scrip for a specific unstructured grid.
Definition: wmscrpmd.F90:551
w3triamd::get_boundary
subroutine get_boundary(MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV, NEIGHBOR_NEXT)
Find boundary points.
Definition: w3triamd.F90:2481
scrip_interface::scrip_clear
subroutine scrip_clear
Definition: scrip_interface.F90:82
wmscrpmd::get_unstructured_vertex_degree
subroutine get_unstructured_vertex_degree(MNP, MNE, TRIGP, TRIGINCD)
This function returns the list of incidences.
Definition: wmscrpmd.F90:1448
scrip_grids::grid2_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lon
Definition: scrip_grids.f:120
scrip_grids::grid1_frac
real(scrip_r8), dimension(:), allocatable, target, save grid1_frac
Definition: scrip_grids.f:103
scrip_remap_vars
Definition: scrip_remap_vars.f:40
wmscrpmd::get_scrip_info
subroutine get_scrip_info(ID_GRD, GRID_CENTER_LON, GRID_CENTER_LAT, GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
Compute grid for scrip for a specific structured grid.
Definition: wmscrpmd.F90:1165
scrip_grids::grid1_units
character(scrip_charlength), save grid1_units
Definition: scrip_grids.f:80
scrip_remap_vars::num_wts
integer(scrip_i4), save num_wts
Definition: scrip_remap_vars.f:66
scrip_grids::grid2_frac
real(scrip_r8), dimension(:), allocatable, target, save grid2_frac
Definition: scrip_grids.f:103
scrip_grids::grid2_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lon
Definition: scrip_grids.f:103