#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE WMSCRPMD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III |
!/ | E. Rogers, M. Dutour, |
!/ | A. Roland, F. Ardhuin |
!/ | FORTRAN 90 |
!/ | Last update : 10-Dec-2014 |
!/ +-----------------------------------+
!/
!/ 06-Sep-2012 : Origination, transfer from WMGRIDMD ( version 4.08 )
!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
!/
!/ Copyright 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 :
!
! Routines to determine and process grid dependencies in the
! multi-grid wave model.
!
! 2. Variables and types :
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! scrip_wrapper Subr. Public as the name says ...
! get_scrip_info_structured Subr. Public as the name says ...
! get_scrip_info_unstructured Subr. Public as the name says ...
! get_scrip_info Subr. Public as the name says ...
! scrip_info_renormalization Subr. Public as the name says ...
! TRIANG_INDEXES Subr. Public as the name says ...
! get_unstructured_vertex_degree Subr. Public as the name says ...
! GET_BOUNDARY Subr. Public as the name says ...
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! get_unstructured_vertex_degree Subr. W3TRIAMD Manage data structures
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! - The subroutines TRIANG_INDEXES, get_unstructured_vertex_degree, and
! GET_BOUNDARY should probably be renamed and moved to the module w3triamd
!
! 6. Switches :
!
!
! !/S Enable subroutine tracing.
! !/Tn Enable test output.
!
! 7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
!/ Specify default accessibility
!/
PUBLIC
!/
!/ Module private variable for checking error returns
!/
INTEGER, PRIVATE :: ISTAT
!/
CONTAINS
!/ ------------------------------------------------------------------- /
SUBROUTINE SCRIP_WRAPPER (ID_SRC, ID_DST, &
MAPSTA_SRC,MAPST2_SRC,FLAGLL,GRIDSHIFT,L_MASTER, &
L_READ,L_TEST)
!/ +-----------------------------------+
!/ | WAVEWATCH III |
!/ | E. Rogers, M. Dutour, A. Roland |
!/ | FORTRAN 90 |
!/ | Last update : 10-Dec-2014 !
!/ +-----------------------------------+
!/
! 1. Original author :
!
! Erick Rogers, NRL
!
! 2. Last update :
!
! See revisions.
!
! 3. Revisions :
!
! 29-Apr-2011 : Origination ( version 4.01 )
! 20-Feb-2012 : Mathieu Dutour Sikiric, Aron Roland Z&P
! Modification is to split the code into several
! subroutines
! ---get_scrip_info_structured (the structured case)
! ---get_scrip_info (chooses according to FD/FE)
! ---scrip_info_renormalization (conv_x and all that)
! 11-Apr-2013 Kevin Lind
! Added code for reading/writing SCRIP remap files
!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
!
! 4. Copyright :
!
! 5. Purpose :
!
! Compute grid information required by SCRIP
!
! 6. Method :
!
! 7. Parameters, Variables and types :
!
! 8. Called by :
!
! Subroutine WMGHGH
!
! 9. Subroutines and functions used :
!
! Subroutine SCRIP
!
! 10. Error messages:
!
! 11. Remarks :
!
! 12. Structure :
!
! 13. Switches :
!
! 14. Source code :
USE SCRIP_GRIDS
USE SCRIP_REMAP_VARS
USE SCRIP_CONSTANTS
USE SCRIP_KINDSMOD
USE SCRIP_INTERFACE
USE W3SERVMD, ONLY: EXTCDE
! USE W3GDATMD, ONLY : GRIDS
IMPLICIT NONE
INTEGER(SCRIP_I4), INTENT(IN) :: ID_SRC, ID_DST
INTEGER(SCRIP_I4), INTENT(IN) :: MAPSTA_SRC(:,:)
INTEGER(SCRIP_I4), INTENT(IN) :: MAPST2_SRC(:,:)
LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: FLAGLL
REAL (SCRIP_R8), INTENT(IN) :: GRIDSHIFT
LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_MASTER ! Am I the master processor (do I/O)?
LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_READ ! Do I read the remap file?
LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_TEST ! Whether to include test output
! in subroutines
!/ ------------------------------------------------------------------- /
!/ local variables
!/
INTEGER(SCRIP_I4) :: IREC,I,J,NI,NJ,IDUM,NK,K, &
ILINK,IW,ICORNER,NGOODPNTS, &
NBADPNTS
INTEGER(SCRIP_I4) :: ISRC,JSRC,KSRC,IPNT,KDST, &
NI_SRC
REAL (SCRIP_R8) :: LAT_CONVERSION,OFFSET
REAL (SCRIP_R8) :: CONV_DX,CONV_DY,WEIGHT
REAL (SCRIP_R8) :: WTSUM
!/T38 CHARACTER (LEN=10) :: CDATE_TIME(3)
!/T38 INTEGER :: DATE_TIME(8)
!/T38 INTEGER :: ELAPSED_TIME, BEG_TIME, &
!/T38 END_TIME
! test output for input variables
!/T38 if(l_master)write(*,*)'flagll = ',flagll
!/T38 if(l_master)write(*,*)'gridshift = ',gridshift
!
! START: universal settings
!
! Set variables for converting to degrees
! notes: SCRIP only operates on spherical coordinates, so for the case
! where the problem is specified by the user as in a
! meters/cartesian coordinate system, it is necessary to make
! a temporary conversion to a "fake" spherical coordinate grid,
! to keep SCRIP happy. The good news here is that multi-grid
! meters-grid simulations will be very rare: we will probably only
! encounter them in the context of simple test cases. Strictly
! speaking, this conversion does not even need to be physically
! correct, e.g. we could say that 1000 km is 1 deg....as long as
! we are consistent between grids.
! Potential future improvement: make conv_dy and offset such that dst grid
! covers a specific longitude range, e.g. 1 deg east to 2 deg east
!
! START: set up src grid
!
!notes: when we work out how to interface with an unstructured grid,
! we will need to revisit this issue of how to set grid1_rank, etc.
! strategy: declare variables in grid module, but allocate them here.
!
GRID1_UNITS='degrees' ! the other option is radians...we don't use this
GRID1_NAME='src' ! this is not used, except for netcdf output
CALL GET_SCRIP_INFO(ID_SRC, &
& GRID1_CENTER_LON, GRID1_CENTER_LAT, &
& GRID1_CORNER_LON, GRID1_CORNER_LAT, GRID1_MASK, &
& GRID1_DIMS, GRID1_SIZE, GRID1_CORNERS, GRID1_RANK)
GRID2_UNITS='degrees'
GRID2_NAME='dst'
CALL GET_SCRIP_INFO(ID_DST, &
& GRID2_CENTER_LON, GRID2_CENTER_LAT, &
& GRID2_CORNER_LON, GRID2_CORNER_LAT, GRID2_MASK, &
& GRID2_DIMS, GRID2_SIZE, GRID2_CORNERS, GRID2_RANK)
IF(FLAGLL)THEN
CONV_DX=ONE
CONV_DY=ONE
OFFSET=ZERO
ELSE
LAT_CONVERSION=ZERO ! lat_conversion
! notes: this is the latitude used for conversion everywhere
! in the grid (approximation)
! (in radians)
! conv_dy=92.6*1200.0 ! physical, =92.6/(3/3600)=111000 m = 111 km
CONV_DY=1.0E+6_SCRIP_R8 ! non-physical, 1e+6=1 deg
CONV_DX=COS(LAT_CONVERSION)*CONV_DY
! notes: offset (in meters), is necessary so that our grid does
! not lie on the branch cut
OFFSET=75000.0_SCRIP_R8-MIN(MINVAL(GRID1_CENTER_LON), &
& MINVAL(GRID2_CENTER_LON))
ENDIF
!.....test output
!/T38 write(*,*)'l_master = ',l_master
!/T38 if(l_master)then
!/T38 write(*,*)'conv_dx=', conv_dx
!/T38 write(*,*)'conv_dy=', conv_dy
!/T38 write(*,*)'offset = ',offset
!/T38 write(*,*)'grid1_size=', grid1_size
!/T38 write(*,*)'grid2_size=', grid2_size
!/T38 write(*,*)'l_read = ',l_read
!/T38 write(*,*)'minval(grid1_center_lon) = ',minval(grid1_center_lon)
!/T38 write(*,*)'maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
!/T38 write(*,*)'minval(grid1_center_lat) = ',minval(grid1_center_lat)
!/T38 write(*,*)'maxval(grid1_center_lat) = ',maxval(grid1_center_lat)
!/T38 write(*,*)'minval(grid2_center_lon) = ',minval(grid2_center_lon)
!/T38 write(*,*)'maxval(grid2_center_lon) = ',maxval(grid2_center_lon)
!/T38 write(*,*)'minval(grid2_center_lat) = ',minval(grid2_center_lat)
!/T38 write(*,*)'maxval(grid2_center_lat) = ',maxval(grid2_center_lat)
!/T38 endif
CALL SCRIP_INFO_RENORMALIZATION( &
& GRID1_CENTER_LON, GRID1_CENTER_LAT, &
& GRID1_CORNER_LON, GRID1_CORNER_LAT, GRID1_MASK, &
& GRID1_DIMS, GRID1_SIZE, GRID1_CORNERS, GRID1_RANK, &
& CONV_DX, CONV_DY, OFFSET, GRIDSHIFT)
CALL SCRIP_INFO_RENORMALIZATION( &
& GRID2_CENTER_LON, GRID2_CENTER_LAT, &
& GRID2_CORNER_LON, GRID2_CORNER_LAT, GRID2_MASK, &
& GRID2_DIMS, GRID2_SIZE, GRID2_CORNERS, GRID2_RANK, &
& CONV_DX, CONV_DY, OFFSET, ZERO)
!.....Set constants for thresholding weights:
FRAC_LOWEST =1.E-3_SCRIP_R8
FRAC_HIGHEST=ONE+1.E-3_SCRIP_R8
WT_LOWEST =ZERO
WT_HIGHEST =ONE+1.E-7_SCRIP_R8
!/T38 call date_and_time (CDATE_TIME(1), CDATE_TIME(2), CDATE_TIME(3), DATE_TIME)
!/T38 beg_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
CALL SCRIP(ID_SRC, ID_DST, L_MASTER, L_READ, L_TEST)
!/T38 call date_and_time (CDATE_TIME(1), CDATE_TIME(2), CDATE_TIME(3), DATE_TIME)
!/T38 end_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
!/T38 elapsed_time = end_time - beg_time
!/T38 write(0,*) "SCRIP: ", elapsed_time, " MSEC"
!/T38 if(l_master)write(*,*)'new minval(grid1_center_lon) = ',minval(grid1_center_lon)
!/T38 if(l_master)write(*,*)'new maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
!.....notes: at this point we have the following useful variables:
! num_wts, e.g. num_wts=3....for first order conservative remapping,
! only the first one is relevant, the other two are for second-order
! remapping.
! max_links_map1, e.g. max_links_map1=1369,
! grid2_size, e.g. grid2_size=1849,
! wts_map1(num_wts,max_links_map1), e.g. wts_map1(3,1369),
! grid1_add_map1(max_links_map1), e.g. grid1_add_map1(1369),
! grid2_add_map1(max_links_map1), e.g. grid2_add_map1(1369),
! grid2_frac(grid2_size), e.g. grid2_frac(1849),
! (see earlier versions for notes re: equivalency in netcdf/matlab)
!
!.....test output (optional)
!
!.....note re: notation: I use k for the combined i/j array, similar to isea,
! but not necessarily the same as isea since some points may
! be land etc.
!/T38 if(l_master)then
!/T38 do k=1,grid2_size
!/T38 write(403,*)grid2_frac(k)
!/T38 end do
!/T38 do ilink=1,max_links_map1
!/T38 write(405,'(999(1x,f20.7))')(wts_map1(iw,ilink),iw=1,num_wts)
!/T38 end do
!/T38 do ilink=1,max_links_map1
!/T38 write(406,'(i20)')grid1_add_map1(ilink) ! equivalent to
!/T38 ! my "src_address"
!/T38 write(407,'(i20)')grid2_add_map1(ilink) ! equivalent to
!/T38 ! my "dst_address"
!/T38 end do
!/T38 endif
!.....organize results and return to wmghgh.
! what we need, for purpose of feeding back to ww3, for each dst grid node:
! a) the set of src grid nodes, in terms of isea, for which
! weights are available
! b) the corresponding set of weights
ALLOCATE ( WGTDATA(GRID2_SIZE), STAT=ISTAT )
CHECK_ALLOC_STATUS ( ISTAT )
!.....step 1: count up NR0, NR1, NR2, NRL, NLOC (NR1 and NLOC are denoted
! "n" here)
! It is especially important to determine how large npnts gets,
! so that we can allocate arrays properly
WGTDATA%NR0=0
WGTDATA%NR2=0
WGTDATA%NRL=0
WGTDATA%N=0
NI_SRC=GRID1_DIMS(1)
NGOODPNTS=0
NBADPNTS=0
DO ILINK=1,MAX_LINKS_MAP1
!........note that this pair of if-thens *must* be consistent with the
!........single if-then below
IF ((GRID2_FRAC(GRID2_ADD_MAP1(ILINK))>FRAC_LOWEST) .AND. &
(GRID2_FRAC(GRID2_ADD_MAP1(ILINK))<FRAC_HIGHEST)) THEN
! then consider this link either to include, or to print a warning
IF( (WTS_MAP1(1,ILINK)>=WT_LOWEST) .AND. &
(WTS_MAP1(1,ILINK)<=WT_HIGHEST) ) THEN
KSRC=GRID1_ADD_MAP1(ILINK)
JSRC=INT((KSRC-1)/NI_SRC)+1
ISRC=KSRC-(JSRC-1)*NI_SRC
IF (MAPSTA_SRC(JSRC,ISRC).EQ.0) THEN ! excluded point
WGTDATA(GRID2_ADD_MAP1(ILINK))%NR0 = &
WGTDATA(GRID2_ADD_MAP1(ILINK))%NR0 + 1
IF (MAPST2_SRC(JSRC,ISRC).EQ.0)THEN
WGTDATA(GRID2_ADD_MAP1(ILINK))%NRL = &
WGTDATA(GRID2_ADD_MAP1(ILINK))%NRL + 1
ENDIF
ELSE IF (ABS(MAPSTA_SRC(JSRC,ISRC)).EQ.1) THEN
! sea point
WGTDATA(GRID2_ADD_MAP1(ILINK))%N= &
WGTDATA(GRID2_ADD_MAP1(ILINK))%N+1
ELSE IF (ABS(MAPSTA_SRC(JSRC,ISRC)).EQ.2) THEN
! bnd point
WGTDATA(GRID2_ADD_MAP1(ILINK))%NR2 = &
WGTDATA(GRID2_ADD_MAP1(ILINK))%NR2 + 1
END IF
NGOODPNTS=NGOODPNTS+1
ELSEIF ( (GRID1_FRAC(GRID1_ADD_MAP1(ILINK))>FRAC_LOWEST) &
.AND. (GRID1_FRAC(GRID1_ADD_MAP1(ILINK))<FRAC_HIGHEST)&
) THEN
! note that we don't bother giving a warning for the
! cases where grid1_frac is strange.
NBADPNTS=NBADPNTS+1
! we also don't bother giving a warning if we've already
! given a lot of warnings (we keep counting, though)
IF((NBADPNTS.LT.5).AND.L_MASTER)THEN ! print warnings
WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
WRITE(*,'(A)')'WARNING: SCRIP weight problem '
WRITE(*,'(4x,A,I7,A,I7)')'ilink = ',ilink,' out of ',&
max_links_map1
WRITE(*,'(4x,A,I7)')'grid1_add_map1(ilink) = ', &
grid1_add_map1(ilink)
WRITE(*,'(4x,A,I7)')'grid2_add_map1(ilink) = ', &
grid2_add_map1(ilink)
WRITE(*,'(4x,A,E12.4)')'wts_map1(1,ilink) = ', &
wts_map1(1,ilink)
WRITE(*,'(4x,A,F10.5)') &
'grid1_frac(grid1_add_map1(ilink)) = ', &
grid1_frac(grid1_add_map1(ilink))
WRITE(*,'(4x,A,F10.5)') &
'grid2_frac(grid2_add_map1(ilink)) = ', &
grid2_frac(grid2_add_map1(ilink))
WRITE(*,'(4x,A,F10.5)')'grid1_center_lat = ', &
grid1_center_lat(grid1_add_map1(ilink))
WRITE(*,'(4x,A,F10.5)')'grid1_center_lon = ', &
grid1_center_lon(grid1_add_map1(ilink))
WRITE(*,'(4x,A,F10.5)')'grid2_center_lat = ', &
grid2_center_lat(grid2_add_map1(ilink))
WRITE(*,'(4x,A,F10.5)')'grid2_center_lon = ', &
grid2_center_lon(grid2_add_map1(ilink))
WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
ENDIF
ENDIF
ENDIF
END DO
IF((NBADPNTS.GT.0).AND.L_MASTER)THEN
WRITE(*,'(4x,A,I5,A)')'We had problems in ',NBADPNTS, &
' points.'
WRITE(*,'(4x,I8,A)')NGOODPNTS,' points appear to be OK.'
ENDIF
IF( (NBADPNTS.GT.(NGOODPNTS/30)) .AND.L_MASTER )THEN
WRITE(*,'(4x,A)')'Error: excessive SCRIP failure. Stopping.'
STOP 'wmscrpmd, case 1'
ENDIF
!.....step 2: allocate according to the size "n" determined above
DO KDST=1,GRID2_SIZE
ALLOCATE ( WGTDATA(KDST)%W(WGTDATA(KDST)%N), STAT=ISTAT )
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE ( WGTDATA(KDST)%K(WGTDATA(KDST)%N), STAT=ISTAT )
CHECK_ALLOC_STATUS ( ISTAT )
WGTDATA(KDST)%N=0
END DO
!.....step 3: save weights
DO ILINK=1,MAX_LINKS_MAP1
KSRC=GRID1_ADD_MAP1(ILINK)
JSRC=INT((KSRC-1)/NI_SRC)+1
ISRC=KSRC-(JSRC-1)*NI_SRC
!........note that this single if-then *must* be consistent with the
!........pair of if-thens above
IF ((GRID2_FRAC(GRID2_ADD_MAP1(ILINK))>FRAC_LOWEST) .AND. &
& (GRID2_FRAC(GRID2_ADD_MAP1(ILINK))<FRAC_HIGHEST) .AND. &
& (WTS_MAP1(1,ILINK)>=WT_LOWEST) .AND. &
& (WTS_MAP1(1,ILINK)<=WT_HIGHEST))THEN
IF (ABS(MAPSTA_SRC(JSRC,ISRC)).EQ.1) THEN ! sea point
WGTDATA(GRID2_ADD_MAP1(ILINK))%N= &
& WGTDATA(GRID2_ADD_MAP1(ILINK))%N+1
WGTDATA(GRID2_ADD_MAP1(ILINK))%W(WGTDATA( &
& GRID2_ADD_MAP1(ILINK))%N)=WTS_MAP1(1,ILINK)
WGTDATA(GRID2_ADD_MAP1(ILINK))%K(WGTDATA( &
& GRID2_ADD_MAP1(ILINK))%N)=GRID1_ADD_MAP1(ILINK)
ENDIF
ENDIF
END DO
!.....step 4: re-normalize weights. This is necessary because we called
!.....scrip without the mask. Now that we have the mask in place, we need
!.....to re-normalize the weights.
DO KDST=1,GRID2_SIZE
IF (WGTDATA(KDST)%N > 0) THEN
WTSUM=ZERO
DO IPNT=1,WGTDATA(KDST)%N
WTSUM=WTSUM+WGTDATA(KDST)%W(IPNT)
END DO
DO IPNT=1,WGTDATA(KDST)%N
WGTDATA(KDST)%W(IPNT)=WGTDATA(KDST)%W(IPNT)/WTSUM
END DO
END IF
END DO
CALL SCRIP_CLEAR
END SUBROUTINE SCRIP_WRAPPER
!/ ------------------------------------------------------------------- /
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)
!/ +-----------------------------------+
!/ | WAVEWATCH III |
!/ | M. Dutour, A. Roland |
!/ | FORTRAN 90 |
!/ | Last update : 10-Dec-2014 !
!/ +-----------------------------------+
!/
! 1. Original author :
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
! 2. Last update :
!
! See revisions.
!
! 3. Revisions :
!
! 20-Feb-2012 : Origination
!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
!
! 4. Copyright :
!
! 5. Purpose :
!
! 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.
!
! 6. Method :
!
! 7. Parameters, Variables and types :
!
! 8. Called by :
!
! Subroutine get_scrip_info
!
! 9. Subroutines and functions used :
!
! 10. Error messages:
!
! 11. Remarks :
!
! 12. Structure :
!
! 13. Switches :
!
! 14. Source code :
USE W3SERVMD, ONLY: EXTCDE
USE W3GDATMD, ONLY : GRIDS
IMPLICIT NONE
INTEGER, INTENT(IN) :: ID_GRD
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LON(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LAT(:)
LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LON(:,:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LAT(:,:)
INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
INTEGER DIRAPPROACH, DUALAPPROACH, THEAPPROACH
INTEGER MNE, MNP, IE, IP, I
INTEGER NBPLUS, NBMINUS
INTEGER I1, I2, I3
REAL*8 :: ELON1, ELON2, ELON3, ELON, ELONC
REAL*8 :: ELAT1, ELAT2, ELAT3, ELAT, ELATC
REAL *8 :: DELTALON12, DELTALON13, DELTALAT12, DELTALAT13
REAL *8 :: THEDET
INTEGER, POINTER :: IOBP(:), TRIGINCD(:)
INTEGER, POINTER :: NEIGHBOR_PREV(:), NEIGHBOR_NEXT(:)
INTEGER, POINTER :: NBASSIGNEDCORNER(:), LISTNBCORNER(:)
INTEGER, POINTER :: STATUS(:), NEXTVERT(:), PREVVERT(:), FINALVERT(:)
INTEGER :: MAXCORNER, NBCORNER
INTEGER :: IDX, IPNEXT, IPPREV, NB, INEXT, IPREV
REAL*8, POINTER :: LON_CENT_TRIG(:), LAT_CENT_TRIG(:)
REAL*8 :: ELONIP, ELONNEXT, ELONPREV, ELONN, ELONP
REAL*8 :: ELATIP, ELATNEXT, ELATPREV, ELATN, ELATP
INTEGER :: ISFINISHED, ZPREV
INTEGER :: DODEBUG
GRID_RANK=2
DIRAPPROACH=1
DUALAPPROACH=2
THEAPPROACH=DUALAPPROACH
MNE=GRIDS(ID_GRD)%NTRI
MNP=GRIDS(ID_GRD)%NX
IF (THEAPPROACH .EQ. DIRAPPROACH) THEN
ALLOCATE(GRID_CENTER_LON(MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CENTER_LAT(MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LON(3,MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LAT(3,MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_MASK(MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
DO IE=1,MNE
I1=GRIDS(ID_GRD)%TRIGP(IE,1)
I2=GRIDS(ID_GRD)%TRIGP(IE,2)
I3=GRIDS(ID_GRD)%TRIGP(IE,3)
ELON1=GRIDS(ID_GRD)%XYB(I1,1)
ELON2=GRIDS(ID_GRD)%XYB(I2,1)
ELON3=GRIDS(ID_GRD)%XYB(I3,1)
ELAT1=GRIDS(ID_GRD)%XYB(I1,2)
ELAT2=GRIDS(ID_GRD)%XYB(I2,2)
ELAT3=GRIDS(ID_GRD)%XYB(I3,2)
ELON=(ELON1 + ELON2 + ELON3)/3
ELAT=(ELAT1 + ELAT2 + ELAT3)/3
GRID_CENTER_LON(IE)=ELON
GRID_CENTER_LAT(IE)=ELAT
GRID_CORNER_LON(1,IE)=ELON1
GRID_CORNER_LON(2,IE)=ELON2
GRID_CORNER_LON(3,IE)=ELON3
GRID_CORNER_LAT(1,IE)=ELAT1
GRID_CORNER_LAT(2,IE)=ELAT2
GRID_CORNER_LAT(3,IE)=ELAT3
GRID_MASK(IE)=.TRUE.
END DO
GRID_CORNERS=3
END IF
IF (THEAPPROACH .EQ. DUALAPPROACH) THEN
ALLOCATE(TRIGINCD(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(IOBP(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(NEIGHBOR_NEXT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(NEIGHBOR_PREV(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(NBASSIGNEDCORNER(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(LISTNBCORNER(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(STATUS(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(NEXTVERT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(PREVVERT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(FINALVERT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(LON_CENT_TRIG(MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(LAT_CENT_TRIG(MNE), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
CALL GET_UNSTRUCTURED_VERTEX_DEGREE (MNP, MNE, &
GRIDS(ID_GRD)%TRIGP, TRIGINCD)
CALL GET_BOUNDARY(MNP, MNE, GRIDS(id_grd)%TRIGP, IOBP, &
NEIGHBOR_PREV, NEIGHBOR_NEXT)
MAXCORNER=0
DO IP=1,MNP
IF (NEIGHBOR_NEXT(IP) .EQ. 0) THEN
NBCORNER=TRIGINCD(IP)
ELSE
NBCORNER=TRIGINCD(IP) + 3
END IF
LISTNBCORNER(IP)=NBCORNER
IF (NBCORNER .GT. MAXCORNER) THEN
MAXCORNER=NBCORNER
END IF
END DO
GRID_CORNERS=MAXCORNER
ALLOCATE(GRID_CENTER_LON(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CENTER_LAT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LON(MAXCORNER,MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LAT(MAXCORNER,MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_MASK(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
NBASSIGNEDCORNER(:)=0
DO IP=1,MNP
GRID_MASK(IP)=.TRUE.
IF (NEIGHBOR_NEXT(IP) .GT. 0) THEN
IPNEXT=NEIGHBOR_NEXT(IP)
IPPREV=NEIGHBOR_PREV(IP)
ELONIP=DBLE(GRIDS(ID_GRD)%XYB(IP,1))
ELATIP=DBLE(GRIDS(ID_GRD)%XYB(IP,2))
ELONNEXT=DBLE(GRIDS(ID_GRD)%XYB(IPNEXT,1))
ELATNEXT=DBLE(GRIDS(ID_GRD)%XYB(IPNEXT,2))
ELONPREV=DBLE(GRIDS(ID_GRD)%XYB(IPPREV,1))
ELATPREV=DBLE(GRIDS(ID_GRD)%XYB(IPPREV,2))
ELONN=(ELONIP+ELONNEXT)/2
ELATN=(ELATIP+ELATNEXT)/2
ELONP=(ELONIP+ELONPREV)/2
ELATP=(ELATIP+ELATPREV)/2
GRID_CORNER_LON(1,IP)=ELONN
GRID_CORNER_LAT(1,IP)=ELATN
GRID_CORNER_LON(2,IP)=ELONIP
GRID_CORNER_LAT(2,IP)=ELATIP
GRID_CORNER_LON(3,IP)=ELONP
GRID_CORNER_LAT(3,IP)=ELATP
NBASSIGNEDCORNER(IP)=3
END IF
END DO
DO IP=1,MNP
GRID_CENTER_LON(IP)=DBLE(GRIDS(ID_GRD)%XYB(IP,1))
GRID_CENTER_LAT(IP)=DBLE(GRIDS(ID_GRD)%XYB(IP,2))
END DO
NBPLUS=0
NBMINUS=0
DO IE=1,MNE
I1=GRIDS(ID_GRD)%TRIGP(IE,1)
I2=GRIDS(ID_GRD)%TRIGP(IE,2)
I3=GRIDS(ID_GRD)%TRIGP(IE,3)
ELON1=DBLE(GRIDS(ID_GRD)%XYB(I1,1))
ELON2=DBLE(GRIDS(ID_GRD)%XYB(I2,1))
ELON3=DBLE(GRIDS(ID_GRD)%XYB(I3,1))
ELAT1=DBLE(GRIDS(ID_GRD)%XYB(I1,2))
ELAT2=DBLE(GRIDS(ID_GRD)%XYB(I2,2))
ELAT3=DBLE(GRIDS(ID_GRD)%XYB(I3,2))
DELTALON12=ELON2 - ELON1
DELTALON13=ELON3 - ELON1
DELTALAT12=ELAT2 - ELAT1
DELTALAT13=ELAT3 - ELAT1
THEDET=DELTALON12*DELTALAT13 - DELTALON13*DELTALAT12
IF (THEDET.GT.0) THEN
NBPLUS=NBPLUS+1
END IF
IF (THEDET.LT.0) THEN
NBMINUS=NBMINUS+1
END IF
ELON=(ELON1 + ELON2 + ELON3)/3
ELAT=(ELAT1 + ELAT2 + ELAT3)/3
LON_CENT_TRIG(IE)=ELON
LAT_CENT_TRIG(IE)=ELAT
END DO
DODEBUG=0
IF (DODEBUG.EQ.1) THEN
print *, 'nbplus=', nbplus, ' nbminus=', nbminus
END IF
STATUS(:) = 0
NEXTVERT(:) = 0
PREVVERT(:) = 0
DO IE=1,MNE
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=GRIDS(ID_GRD)%TRIGP(IE,I)
IPNEXT=GRIDS(ID_GRD)%TRIGP(IE,INEXT)
IPPREV=GRIDS(ID_GRD)%TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
IF (NEIGHBOR_NEXT(IP).EQ.0) THEN
STATUS(IP)=1
FINALVERT(IP)=IPPREV
PREVVERT(IP)=IPPREV
NEXTVERT(IP)=IPNEXT
ELSE
IF (NEIGHBOR_PREV(IP).EQ.IPPREV) THEN
STATUS(IP)=1
PREVVERT(IP)=IPPREV
NEXTVERT(IP)=IPNEXT
FINALVERT(IP)=NEIGHBOR_NEXT(IP)
END IF
END IF
END IF
END DO
END DO
STATUS(:)=0
DO
ISFINISHED=1
DO IE=1,MNE
ELON=LON_CENT_TRIG(IE)
ELAT=LAT_CENT_TRIG(IE)
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=GRIDS(ID_GRD)%TRIGP(IE,I)
IPNEXT=GRIDS(ID_GRD)%TRIGP(IE,INEXT)
IPPREV=GRIDS(ID_GRD)%TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
ISFINISHED=0
ZPREV=PREVVERT(IP)
IF (ZPREV.EQ.IPPREV) THEN
IDX=NBASSIGNEDCORNER(IP)
IDX=IDX+1
GRID_CORNER_LON(IDX,IP)=ELON
GRID_CORNER_LAT(IDX,IP)=ELAT
NBASSIGNEDCORNER(IP)=IDX
PREVVERT(IP)=IPNEXT
IF (IPNEXT.EQ.FINALVERT(IP)) THEN
STATUS(IP)=1
END IF
END IF
END IF
END DO
END DO
IF (ISFINISHED.EQ.1) THEN
EXIT
END IF
END DO
DO IP=1,MNP
IF (NBASSIGNEDCORNER(IP).NE.LISTNBCORNER(IP)) THEN
WRITE(*,*) 'Incoherent number at IP=', IP
WRITE(*,*) ' NbAssignedCorner(IP)=', NbAssignedCorner(IP)
WRITE(*,*) ' ListNbCorner(IP)=', ListNbCorner(IP)
WRITE(*,*) ' N_N=', NEIGHBOR_NEXT(IP), 'N_P=', NEIGHBOR_PREV(IP)
WRITE(*,*) ' TrigIncd=', TrigIncd(IP)
STOP 'wmscrpmd, case 2'
END IF
END DO
! if the number of corner is below threshold, we have to
! add some more.
DO IP=1,MNP
NB=NBASSIGNEDCORNER(IP)
IF (NB .LT. MAXCORNER) THEN
ELON=GRID_CORNER_LON(NB,IP)
ELAT=GRID_CORNER_LAT(NB,IP)
DO IDX=NB+1,MAXCORNER
GRID_CORNER_LON(IDX,IP)=ELON
GRID_CORNER_LAT(IDX,IP)=ELAT
END DO
END IF
END DO
DEALLOCATE(NBASSIGNEDCORNER, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(LISTNBCORNER, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(TRIGINCD, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(IOBP, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(NEIGHBOR_PREV, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(NEIGHBOR_NEXT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(STATUS, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(NEXTVERT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(PREVVERT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(FINALVERT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(LON_CENT_TRIG, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(LAT_CENT_TRIG, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_DIMS(2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
GRID_DIMS(1)=MNP
GRID_DIMS(2)=1
GRID_SIZE=MNP
END IF
END SUBROUTINE GET_SCRIP_INFO_UNSTRUCTURED
!/ ------------------------------------------------------------------- /
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)
!/ +-----------------------------------+
!/ | WAVEWATCH III |
!/ | M. Dutour, A. Roland |
!/ | FORTRAN 90 |
!/ | Last update : 10-Dec-2014 !
!/ +-----------------------------------+
!/
! 1. Original author :
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
! 2. Last update :
!
! See revisions.
!
! 3. Revisions :
!
! 20-Feb-2012 : Origination, this is adapted from Erick Rogers
! code by splitting the code into sections.
!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
!
! 4. Copyright :
!
! 5. Purpose :
!
! 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
!
! 6. Method :
!
! 7. Parameters, Variables and types :
!
! 8. Called by :
!
! Subroutine get_scrip_info
!
! 9. Subroutines and functions used :
!
! 10. Error messages:
!
! 11. Remarks :
!
! 12. Structure :
!
! 13. Switches :
!
! 14. Source code :
USE W3SERVMD, ONLY: EXTCDE
USE W3GDATMD, ONLY : GRIDS
USE SCRIP_CONSTANTS, ONLY : HALF
IMPLICIT NONE
INTEGER, INTENT(IN) :: ID_GRD
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LON(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LAT(:)
LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LON(:,:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LAT(:,:)
INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
!
REAL*8, ALLOCATABLE :: XIN_GRD(:,:), YIN_GRD(:,:)
REAL*8, ALLOCATABLE :: DXDP_GRD(:,:), DXDQ_GRD(:,:)
REAL*8, ALLOCATABLE :: DYDP_GRD(:,:), DYDQ_GRD(:,:)
INTEGER :: N1, N2, NI, NJ
INTEGER :: IREC, J, I
GRID_RANK=2
N1=SIZE(GRIDS(ID_GRD)%XGRD,1)
N2=SIZE(GRIDS(ID_GRD)%XGRD,2)
ALLOCATE(XIN_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(YIN_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(DXDP_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(DXDQ_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(DYDP_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(DYDQ_GRD(N1,N2), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
XIN_GRD=DBLE(GRIDS(ID_GRD)%XGRD)
YIN_GRD=DBLE(GRIDS(ID_GRD)%YGRD)
DXDP_GRD=DBLE(GRIDS(ID_GRD)%DXDP)
DXDQ_GRD=DBLE(GRIDS(ID_GRD)%DXDQ)
DYDP_GRD=DBLE(GRIDS(ID_GRD)%DYDP)
DYDQ_GRD=DBLE(GRIDS(ID_GRD)%DYDQ)
ALLOCATE(GRID_DIMS(GRID_RANK), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
GRID_DIMS(1)=N2
NI=N2
GRID_DIMS(2)=N1
NJ=N1
GRID_SIZE=NI*NJ ! hardwired: logically rectangular grid
GRID_CORNERS=4 ! hardwired: each cell has 4 corners
!.....notes: unfortunately, scrip only works for spherical coordinates.
! thus, if we want to have a multi-grid case in meters, we have to
! fake it. fortunately, it should be pretty rare to have a multi-grid
! case in meters.
ALLOCATE(GRID_CENTER_LON(NI*NJ), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CENTER_LAT(NI*NJ), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LON(4,NI*NJ), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_CORNER_LAT(4,NI*NJ), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(GRID_MASK(NI*NJ), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
!.....notes: this "gridshift" variable is included because SCRIP sometimes
! has trouble when grids cell locations are identical between
! the two grids. Thus we apply this to one of the two grids.
!.....notes: The following block of code could be converted to a subroutine.
! Since it is called twice, this would save a little space.
IREC=0
DO J=1,NJ
DO I=1,NI
IREC=IREC+1
GRID_CENTER_LON(IREC)=XIN_GRD(J,I)
GRID_CENTER_LAT(IREC)=YIN_GRD(J,I)
GRID_MASK(IREC)=.TRUE.
!..notes: normally, we'd apply the mask like this:
! if(abs(mapsta_src(j,i)).eq.1)then
! grid1_mask(irec)=.true.
! else
! grid1_mask(irec)=.false.
! endif
!..but unfortunately, WMGHGH needs information about the overlaying high-res
! cells, even those that are masked, for calculating
! NRL, NR0, NR1, NR2.
!...........corner 1 : halfway to i-1,j-1
GRID_CORNER_LON(1,IREC)=GRID_CENTER_LON(IREC)- &
& HALF*DXDP_GRD(J,I)-HALF*DXDQ_GRD(J,I)
GRID_CORNER_LAT(1,IREC)=GRID_CENTER_LAT(IREC)- &
& HALF*DYDP_GRD(J,I)-HALF*DYDQ_GRD(J,I)
!...........corner 2: halfway to i+1,j-1
GRID_CORNER_LON(2,IREC)=GRID_CENTER_LON(IREC)+ &
& HALF*DXDP_GRD(J,I)-HALF*DXDQ_GRD(J,I)
GRID_CORNER_LAT(2,IREC)=GRID_CENTER_LAT(IREC)+ &
& HALF*DYDP_GRD(J,I)-HALF*DYDQ_GRD(J,I)
!...........corner 3: halfway to i+1,j+1
GRID_CORNER_LON(3,IREC)=GRID_CENTER_LON(IREC)+ &
& HALF*DXDP_GRD(J,I)+HALF*DXDQ_GRD(J,I)
GRID_CORNER_LAT(3,IREC)=GRID_CENTER_LAT(IREC)+ &
& HALF*DYDP_GRD(J,I)+HALF*DYDQ_GRD(J,I)
!...........corner 4: halfway to i-1,j+1
GRID_CORNER_LON(4,IREC)=GRID_CENTER_LON(IREC)- &
& HALF*DXDP_GRD(J,I)+HALF*DXDQ_GRD(J,I)
GRID_CORNER_LAT(4,IREC)=GRID_CENTER_LAT(IREC)- &
& HALF*DYDP_GRD(J,I)+HALF*DYDQ_GRD(J,I)
END DO
END DO
END SUBROUTINE GET_SCRIP_INFO_STRUCTURED
!/ ------------------------------------------------------------------- /
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)
! 1. Original author :
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
! 2. Last update :
!
! See revisions.
!
! 3. Revisions :
!
! 20-Feb-2012 : Origination, this is adapted from Erick Rogers
! code by splitting the code into sections.
!
! 4. Copyright :
!
! 5. Purpose :
!
! Compute grid for scrip for a specific structured grid.
! This is adapted from Erick Rogers code by making it cleaner.
!
! 6. Method :
!
! 7. Parameters, Variables and types :
!
! 8. Called by :
!
! Subroutine scrip_wrapper
!
! 9. Subroutines and functions used :
!
! 10. Error messages:
!
! 11. Remarks :
!
! 12. Structure :
!
! 13. Switches :
!
! 14. Source code :
USE W3SERVMD, ONLY: EXTCDE
USE W3GDATMD, ONLY : GRIDS, UNGTYPE
IMPLICIT NONE
INTEGER, INTENT(IN) :: ID_GRD
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LON(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LAT(:)
LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LON(:,:)
REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LAT(:,:)
INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
REAL*8 :: DLON1, DLAT1, DLON2, DLAT2, THEDET
INTEGER :: I, J
INTEGER :: IC, JC, IP, CHECKSIGNS, NBPLUS, NBMINUS
INTEGER :: PRINTDATA, PRINTMINMAX
REAL*8 :: MINLON, MINLAT, MAXLON, MAXLAT
REAL*8 :: MINLONCORNER, MAXLONCORNER, MINLATCORNER, MAXLATCORNER
IF (GRIDS(ID_GRD)%GTYPE .EQ. UNGTYPE) THEN
CALL 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)
ELSE
CALL 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)
END IF
CHECKSIGNS=0
IF (CHECKSIGNS.EQ.1) THEN
NBPLUS=0
NBMINUS=0
DO IP=1,GRID_SIZE
DO IC=1,GRID_CORNERS
IF (IC.EQ.GRID_CORNERS) THEN
JC=1
ELSE
JC=IC+1
END IF
DLON1=GRID_CORNER_LON(IC,IP)-GRID_CENTER_LON(IP)
DLON2=GRID_CORNER_LON(JC,IP)-GRID_CENTER_LON(IP)
DLAT1=GRID_CORNER_LAT(IC,IP)-GRID_CENTER_LAT(IP)
DLAT2=GRID_CORNER_LAT(JC,IP)-GRID_CENTER_LAT(IP)
THEDET=DLON1*DLAT2 - DLON2*DLAT1
IF (THEDET.GT.0) THEN
NBPLUS=NBPLUS+1
END IF
IF (THEDET.LT.0) THEN
NBMINUS=NBMINUS+1
END IF
END DO
END DO
WRITE(*,*) 'SI nbPlus=', nbPlus, ' nbMinus=', nbMinus
END IF
END SUBROUTINE GET_SCRIP_INFO
!/ ------------------------------------------------------------------- /
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)
! 1. Original author :
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
! Adapted from Erick Rogers scrip_wrapper
!
! 2. Last update :
!
! See revisions.
!
! 3. Revisions :
!
! 20-Feb-2012 : Origination
!
! 4. Copyright :
!
! 5. Purpose :
!
! 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
!
! 6. Method :
!
! We apply various transformations to the longitude latitude
! following here the transformations that were done only in
! finite difference meshes.
!
! 7. Parameters, Variables and types :
!
! 8. Called by :
!
! Subroutine WMGHGH
!
! 9. Subroutines and functions used :
!
! 10. Error messages:
!
! 11. Remarks :
!
! 12. Structure :
!
! 13. Switches :
!
! 14. Source code :
IMPLICIT NONE
REAL*8, INTENT(INOUT) :: GRID_CENTER_LON(:)
REAL*8, INTENT(INOUT) :: GRID_CENTER_LAT(:)
LOGICAL, INTENT(IN) :: GRID_MASK(:)
REAL*8, INTENT(INOUT) :: GRID_CORNER_LON(:,:)
REAL*8, INTENT(INOUT) :: GRID_CORNER_LAT(:,:)
INTEGER, INTENT(IN) :: GRID_DIMS(:)
INTEGER, INTENT(IN) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
REAL*8 :: CONV_DX, CONV_DY, OFFSET, GRIDSHIFT
REAL*8 DEG2RAD
!
INTEGER :: I, J, IP
REAL*8 :: MINLON, MINLAT, MAXLON, MAXLAT, HLON, HLAT
REAL*8 :: MINLONCORNER, MAXLONCORNER, MINLATCORNER, MAXLATCORNER
DO I=1,GRID_SIZE
GRID_CENTER_LON(I)=(GRID_CENTER_LON(I)+OFFSET)/CONV_DX + &
& GRIDSHIFT
GRID_CENTER_LAT(I)=GRID_CENTER_LAT(I)/CONV_DY + &
& GRIDSHIFT
IF(GRID_CENTER_LON(I)>360.0) THEN
GRID_CENTER_LON(I)=GRID_CENTER_LON(I)-360.0
END IF
IF(GRID_CENTER_LON(I)<000.0) THEN
GRID_CENTER_LON(I)=GRID_CENTER_LON(I)+360.0
END IF
DO J=1,GRID_CORNERS
GRID_CORNER_LON(J, I)=(GRID_CORNER_LON(J, I)+OFFSET)/CONV_DX+ &
& GRIDSHIFT
GRID_CORNER_LAT(J, I)=GRID_CORNER_LAT(J, I)/CONV_DY + &
& GRIDSHIFT
IF(GRID_CORNER_LON(J,I)>360.0) THEN
GRID_CORNER_LON(J,I)=GRID_CORNER_LON(J,I)-360.0
END IF
IF(GRID_CORNER_LON(J,I)<000.0) THEN
GRID_CORNER_LON(J,I)=GRID_CORNER_LON(J,I)+360.0
END IF
END DO
END DO
END SUBROUTINE SCRIP_INFO_RENORMALIZATION
!/ ------------------------------------------------------------------- /
SUBROUTINE TRIANG_INDEXES(I, INEXT, IPREV)
! 1. Original author :
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
INTEGER, INTENT(IN) :: I
INTEGER, INTENT(OUT) :: INEXT, IPREV
IF (I.EQ.1) THEN
INEXT=3
ELSE
INEXT=I-1
END IF
IF (I.EQ.3) THEN
IPREV=1
ELSE
IPREV=I+1
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE GET_UNSTRUCTURED_VERTEX_DEGREE (MNP, MNE, TRIGP, &
& TRIGINCD)
! Written:
!
! 20-Feb-2012
!
! Author:
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
! Parameters:
! Input:
! MNP: number of nodes
! INE: list of nodes
! MNE: number of triangles
! Output:
! TrigIncd (number of triangles contained by vertices
!
! Description:
! this function returns the list of incidences
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: MNP, MNE
INTEGER, INTENT(IN) :: TRIGP(:,:)
INTEGER, INTENT(OUT) :: TRIGINCD(:)
INTEGER :: IP, IE, I
TRIGINCD=0
DO IE=1,MNE
DO I=1,3
IP=TRIGP(IE,I)
TRIGINCD(IP)=TRIGINCD(IP) + 1
END DO
END DO
END SUBROUTINE GET_UNSTRUCTURED_VERTEX_DEGREE
!/ ------------------------------------------------------------------- /
SUBROUTINE GET_BOUNDARY(MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV, &
& NEIGHBOR_NEXT)
!/ +-----------------------------------+
!/ | WAVEWATCH III |
!/ | M. Dutour, A. Roland |
!/ | FORTRAN 90 |
!/ | Last update : 10-Dec-2014 !
!/ +-----------------------------------+
!/
! Written:
!
! 20-Feb-2012
!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 )
!
! Author:
!
! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
!
! Parameters:
! Input:
! MNP: number of nodes
! TRIGP: list of nodes
! MNE: number of triangles
! Output:
! NEIGHBOR
!
! Description:
! 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.
!
USE W3SERVMD, ONLY: EXTCDE
IMPLICIT NONE
INTEGER, INTENT(IN) :: MNP, MNE, TRIGP(MNE,3)
INTEGER, INTENT(INOUT) :: IOBP(MNP)
INTEGER, INTENT(INOUT) :: NEIGHBOR_PREV(MNP)
INTEGER, INTENT(INOUT) :: NEIGHBOR_NEXT(MNP)
INTEGER, POINTER :: STATUS(:)
INTEGER, POINTER :: COLLECTED(:)
INTEGER, POINTER :: NEXTVERT(:)
INTEGER, POINTER :: PREVVERT(:)
INTEGER :: IE, I, IP, IP2, IP3
INTEGER :: ISFINISHED, INEXT, IPREV
INTEGER :: IPNEXT, IPPREV, ZNEXT, ZPREV
ALLOCATE(STATUS(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(COLLECTED(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(PREVVERT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
ALLOCATE(NEXTVERT(MNP), STAT=ISTAT)
CHECK_ALLOC_STATUS ( ISTAT )
IOBP = 0
NEIGHBOR_NEXT = 0
NEIGHBOR_PREV = 0
! Now computing the next items
STATUS = 0
NEXTVERT = 0
PREVVERT = 0
DO IE=1,MNE
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=TRIGP(IE,I)
IPNEXT=TRIGP(IE,INEXT)
IPPREV=TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
STATUS(IP)=1
PREVVERT(IP)=IPPREV
NEXTVERT(IP)=IPNEXT
END IF
END DO
END DO
STATUS(:)=0
DO
COLLECTED(:)=0
DO IE=1,MNE
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=TRIGP(IE,I)
IPNEXT=TRIGP(IE,INEXT)
IPPREV=TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
ZNEXT=NEXTVERT(IP)
IF (ZNEXT.EQ.IPPREV) THEN
COLLECTED(IP)=1
NEXTVERT(IP)=IPNEXT
IF (NEXTVERT(IP).EQ.PREVVERT(IP)) THEN
STATUS(IP)=1
END IF
END IF
END IF
END DO
END DO
ISFINISHED=1
DO IP=1,MNP
IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN
STATUS(IP)=-1
NEIGHBOR_NEXT(IP)=NEXTVERT(IP)
END IF
IF (STATUS(IP).EQ.0) THEN
ISFINISHED=0
END IF
END DO
IF (ISFINISHED.EQ.1) THEN
EXIT
END IF
END DO
! Now computing the prev items
STATUS = 0
NEXTVERT = 0
PREVVERT = 0
DO IE=1,MNE
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=TRIGP(IE,I)
IPNEXT=TRIGP(IE,INEXT)
IPPREV=TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
STATUS(IP)=1
PREVVERT(IP)=IPPREV
NEXTVERT(IP)=IPNEXT
END IF
END DO
END DO
STATUS(:)=0
DO
COLLECTED(:)=0
DO IE=1,MNE
DO I=1,3
CALL TRIANG_INDEXES(I, INEXT, IPREV)
IP=TRIGP(IE,I)
IPNEXT=TRIGP(IE,INEXT)
IPPREV=TRIGP(IE,IPREV)
IF (STATUS(IP).EQ.0) THEN
ZPREV=PREVVERT(IP)
IF (ZPREV.EQ.IPNEXT) THEN
COLLECTED(IP)=1
PREVVERT(IP)=IPPREV
IF (PREVVERT(IP).EQ.NEXTVERT(IP)) THEN
STATUS(IP)=1
END IF
END IF
END IF
END DO
END DO
ISFINISHED=1
DO IP=1,MNP
IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN
STATUS(IP)=-1
NEIGHBOR_PREV(IP)=PREVVERT(IP) ! new code
END IF
IF (STATUS(IP).EQ.0) THEN
ISFINISHED=0
END IF
END DO
IF (ISFINISHED.EQ.1) THEN
EXIT
END IF
END DO
! Now making checks
DO IP=1,MNP
IP2=NEIGHBOR_NEXT(IP)
IF (IP2.GT.0) THEN
IP3=NEIGHBOR_PREV(IP2)
IF (ABS(IP3 - IP).GT.0) THEN
WRITE(*,*) 'IP=', IP, ' IP2=', IP2, ' IP3=', IP3
WRITE(*,*) 'We have a dramatic inconsistency'
STOP 'wmscrpmd, case 3'
END IF
END IF
END DO
! Now assigning the boundary IOBP array
DO IP=1,MNP
IF (STATUS(IP).EQ.-1 .AND. IOBP(IP) .EQ. 0) THEN
IOBP(IP)=1
END IF
END DO
DEALLOCATE(STATUS, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(COLLECTED, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(NEXTVERT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
DEALLOCATE(PREVVERT, STAT=ISTAT)
CHECK_DEALLOC_STATUS ( ISTAT )
END SUBROUTINE
!/
!/ End of module WMSCRPMD -------------------------------------------- /
!/
END MODULE WMSCRPMD