#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))=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_LOWEST) .AND. & & (GRID2_FRAC(GRID2_ADD_MAP1(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