Merge remote-tracking branch 'EMC/develop' into NetCDFcomplink

.gitignore

@@ -45,6 +45,10 @@ manual/*.toc
 manual/*.out
 manual/*.dvi
 manual/*.pdf
+model?
+model??
+regtests/list*
+regtests/before
 regtests/matrix*
 regtests/*/work*
 regtests/*/input*/*.nc

model/bin/make_makefile.sh

@@ -887,7 +887,17 @@
                prop=
              source="w3parall w3triamd $stx $nlx $btx $is $uostmd"
                  IO='w3iogrmd'
-                aux="constants w3servmd w3arrymd w3dispmd w3gsrumd w3timemd w3nmlgridmd $pdlibyow $memcode" ;;
+                aux="constants w3servmd w3arrymd w3dispmd w3gsrumd w3timemd w3nmlgridmd $pdlibyow $memcode"
+                if [ "$scrip" = 'SCRIP' ]
+                then
+                  aux="$aux scrip_constants scrip_grids scrip_iounitsmod"
+                  aux="$aux scrip_remap_vars scrip_timers scrip_errormod scrip_interface"
+                  aux="$aux scrip_kindsmod scrip_remap_conservative wmscrpmd" 
+                fi 
+                if [ "$scripnc" = 'SCRIPNC' ]
+                then
+                  aux="$aux scrip_netcdfmod scrip_remap_write scrip_remap_read"
+                fi ;;
      ww3_strt) IDstring='Initial conditions program'
                core=
                data="$memcode w3gdatmd w3wdatmd w3adatmd w3idatmd w3odatmd"

model/bin/w3_make

@@ -135,8 +135,7 @@
 
 # GRID STRT BOUND OUTF OUTP TRCK GRIB GINT GX* UPRSTR LIBWW3
 #  reg_programs = non-NetCDF programs
-  reg_programs="ww3_grid" 
-  reg_programs="$reg_programs ww3_strt"
+  reg_programs="ww3_strt"
   reg_programs="$reg_programs ww3_bound"
   reg_programs="$reg_programs ww3_outf"
   reg_programs="$reg_programs ww3_outp"
@@ -171,6 +170,7 @@
   # MULTI MULTI_ESMF SBS1
   if [ -n "`grep SCRIPNC $switch_file`" ] || [ -n "`grep OASIS $switch_file`" ]  || [ -n "`grep PDLIB $switch_file`" ]
   then
+    cdf_programs="$cdf_programs ww3_grid"
     cdf_programs="$cdf_programs ww3_multi"
     cdf_programs="$cdf_programs ww3_sbs1"
     cdf_programs="$cdf_programs libww3"
@@ -179,6 +179,7 @@
       cdf_programs="$cdf_programs ww3_multi_esmf"
     fi
   else
+    reg_programs="$reg_programs ww3_grid"
     reg_programs="$reg_programs ww3_multi"
     reg_programs="$reg_programs ww3_sbs1"
     reg_programs="$reg_programs libww3"

model/bin/ww3_from_ftp.sh

@@ -8,7 +8,7 @@
 curr_dir=`pwd`
 
 # Set WW3 code version
-ww3ver=v7.07
+ww3ver=v7.12.1
 
 #Get top level directory of ww3 from user: 
 echo -e "\n\n This script will download data from the ftp for WAVEWATCH III "
@@ -51,7 +51,7 @@ cp -r data_regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM regtests/ww
 cp -r data_regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM6
 cp -r data_regtests/ww3_tp2.14/input/toy/*.nc regtests/ww3_tp2.14/input/toy/
 cp -r data_regtests/ww3_tp2.17/input/*     regtests/ww3_tp2.17/input/
-
+cp -r data_regtests/ww3_tp2.21/input/*     regtests/ww3_tp2.21/input/
 
 #Do you want to clean up (aka delete tar file, delete the data_regtests directory) 
 echo -e "\n\n Do you want to delete the tar file ww3_from_ftp.${ww3ver}.tar.gz [y|n]: "

model/ftn/PDLIB/yowpdlibmain.ftn

@@ -1145,6 +1145,7 @@ module yowpdlibMain
     implicit none
     integer I1, I2, I3, stat, IE, NI(3)
     real(rkind) :: DXP1, DXP2, DXP3, DYP1, DYP2, DYP3, DBLTMP, TRIA03
+    logical     :: CROSSES_DATELINE
 
     allocate(PDLIB_SI(npa), PDLIB_CCON(npa), PDLIB_IEN(6,ne), PDLIB_TRIA(ne), stat=stat)
     if(stat/=0) call parallel_abort('SI allocation failure')
@@ -1163,6 +1164,12 @@ module yowpdlibMain
       DYP2=y(I3) - y(I2)
       DXP3=x(I1) - x(I3)
       DYP3=y(I1) - y(I3)
+      CALL ELEMENT_CROSSES_DATELINE(DXP1, DXP2, DXP3, CROSSES_DATELINE)
+      IF (CROSSES_DATELINE) THEN
+        CALL CORRECT_DX_GT180(DXP1)
+        CALL CORRECT_DX_GT180(DXP2)
+        CALL CORRECT_DX_GT180(DXP3)
+      ENDIF
 
       PDLIB_IEN(1,IE) = - DYP2
       PDLIB_IEN(2,IE) =   DXP2
@@ -1185,6 +1192,41 @@ module yowpdlibMain
   end subroutine
 !**********************************************************************
 !*                                                                    *
+!**********************************************************************
+  subroutine ELEMENT_CROSSES_DATELINE(RX1, RX2, RX3, CROSSES_DATELINE)
+!   Purpose: understanding if an element crosses the dateline.
+!   An element crossing the dateline has, e.g. a node with lon < 180
+!   and another 2 with lon > -180
+    IMPLICIT NONE
+    REAL(rkind),  INTENT(IN)  :: RX1, RX2, RX3
+    LOGICAL, INTENT(OUT) :: CROSSES_DATELINE
+    INTEGER :: R1GT180, R2GT180, R3GT180
+    R1GT180 = MERGE(1, 0, ABS(RX1).GT.180)
+    R2GT180 = MERGE(1, 0, ABS(RX2).GT.180)
+    R3GT180 = MERGE(1, 0, ABS(RX3).GT.180)
+! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the dateline
+! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole
+! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline
+    CROSSES_DATELINE = R1GT180+R2GT180+R3GT180 .EQ. 2
+  end subroutine
+!**********************************************************************
+!*                                                                    *
+!**********************************************************************
+  subroutine CORRECT_DX_GT180(DXP)
+!   Purpose: the absolute zonal distance between 2 points is always <= 180
+!            This subroutine corrects the zonal distance to satifsy
+!            this requirement
+    IMPLICIT NONE
+    REAL(rkind), INTENT(INOUT) :: DXP
+    IF (DXP .le. -180) THEN
+      DXP=DXP + 360
+    END IF
+    IF (DXP .ge. 180) THEN
+      DXP=DXP - 360
+    END IF
+  end subroutine
+!**********************************************************************
+!*                                                                    *
 !**********************************************************************
   subroutine ComputeIA_JA_POSI_NNZ
     use yowElementpool, only: ne, ne_global, INE, ielg

model/ftn/SCRIP/scrip_grids.f

@@ -96,6 +96,10 @@ module scrip_grids
      &     special_polar_cell1,       ! cell with only 1 corner at pole
      &     special_polar_cell2        !
 
+      integer (SCRIP_i4), dimension(:), allocatable, target,save ::
+     &             grid1_imask,        ! flag which cells participate
+     &             grid2_imask         ! flag which cells participate
+
       real (SCRIP_r8), dimension(:), allocatable, target, save ::
      &             grid1_center_lat,  ! lat/lon coordinates for
      &             grid1_center_lon,  ! each grid center in radians

model/ftn/w3iopomd.ftn

@@ -242,7 +242,7 @@
       INTEGER                 :: IX1, IY1, IXS, IYS
 !/S      INTEGER, SAVE           :: IENT = 0
 !/O7a      INTEGER                 :: IX0, IXN, IY0, IYN, NNX,         &
-!/O7a                                 KX, KY, JX, IIX
+!/O7a                                 KX, KY, JX, IIX, IX2, IY2, IS1
       INTEGER                 :: IX(4), IY(4)   ! Indices of points used in interp.
       REAL                    :: RD(4)          ! Interpolation coefficient 
 !/O7a      REAL                    :: RD1, RD2, RDTOT, ZBOX(4), DEPTH
@@ -392,7 +392,7 @@
 !/O7a                         RD(4)*ZBOX(4) ) / RDTOT
 !/O7a            WRITE (SCREEN,943) DEPTH
 !
-!/O7a       *** implementation of O7a option with curvilinear grids is incomplete ***
+!/O7a      ! *** implementation of O7a option with curvilinear grids is incomplete ***
 !
 !/O7a            IF ( RD1 .LT. 0.05 ) IX2 = IX1
 !/O7a            IF ( RD1 .GT. 0.95 ) IX1 = IX2

model/ftn/w3triamd.ftn

@@ -614,7 +614,6 @@ CONTAINS
 !
 !  1. purpose: defines open boundary points based on depth
 !  2. Method : a boundary node has more node around it than triangles
-!      
 !    
 !
 !  3. Parameters :
@@ -739,6 +738,8 @@ CONTAINS
          REAL              :: TL1, TL2, TL3, TMPTRIGP
          INTEGER           :: I1, I2, I3
          INTEGER           :: K
+         REAL*8            :: PT(3,2)
+
 !/S        INTEGER                      ::  IENT = 0          
 !/ ------------------------------------------------------------------- /   
 !/S      CALL STRACE (IENT, 'SPATIAL_GRID')   
@@ -748,13 +749,16 @@ CONTAINS
                    I1 = TRIGP(K,1)
                    I2 = TRIGP(K,2)
                    I3 = TRIGP(K,3)
+
+                   CALL FIX_PERIODCITY(I1,I2,I3,XYB,PT) 
 !
 ! cross product of edge-vector  (orientated anticlockwise)
 !                                   
-                   TRIA(K) = REAL( (XYB(I2,2)-XYB(I1,2))      &     !  (Y2-Y1)
-                                  *(XYB(I1,1)-XYB(I3,1))      &     ! *(X1-X3)
-                                  +(XYB(I3,2)-XYB(I1,2))      &     !  (Y3-Y1)*(X2-X1)
-                                  *(XYB(I2,1)-XYB(I1,1))      )*0.5
+
+                   TRIA(K) = REAL( (PT(2,2)-PT(1,2))      &     !  (Y2-Y1)
+                                  *(PT(1,1)-PT(3,1))      &     ! *(X1-X3)
+                                  +(PT(3,2)-PT(1,2))      &     !  (Y3-Y1)*(X2-X1)
+                                  *(PT(2,1)-PT(1,1))      )*0.5
 !
 ! test on negative triangle area, which means that the orientiation is not as assumed to be anticw. 
 ! therefore we swap the nodes !!! 
@@ -832,6 +836,7 @@ CONTAINS
          REAL*8    :: N1(2), N2(2), N3(2)
 	 REAL*8    :: TMP(3)
 	 REAL*8    :: TMPINV(3)
+         REAL*8    :: PT(3,2)
 !/S        INTEGER                      ::  IENT = 0 	 
 !/ ------------------------------------------------------------------- /   
 !/S      CALL STRACE (IENT, 'NVECTRI')   
@@ -845,12 +850,14 @@ CONTAINS
             I2 = TRIGP(IE,2)
             I3 = TRIGP(IE,3)
 
-            P1(1) = XYB(I1,1)
-            P1(2) = XYB(I1,2)
-            P2(1) = XYB(I2,1)
-            P2(2) = XYB(I2,2)
-            P3(1) = XYB(I3,1)
-            P3(2) = XYB(I3,2)
+            CALL FIX_PERIODCITY(I1,I2,I3,XYB,PT) 
+
+            P1(1) = PT(1,1)
+            P1(2) = PT(1,2)
+            P2(1) = PT(2,1)
+            P2(2) = PT(2,2)
+            P3(1) = PT(3,1)
+            P3(2) = PT(3,2)
 !
 ! I1 -> I2, I2 -> I3, I3 -> I1 (anticlockwise orientation is preserved)
 !                    
@@ -1418,6 +1425,7 @@ END SUBROUTINE
      DOUBLE PRECISION             :: x1, x2, x3
      DOUBLE PRECISION             :: y1, y2, y3
      DOUBLE PRECISION             :: s1, s2, s3, sg1, sg2, sg3
+     REAL*8                       :: PT(3,2)
      INTEGER                      :: ITRI
      INTEGER                      :: I1, I2, I3
      INTEGER                      :: nbFound
@@ -1433,15 +1441,17 @@ END SUBROUTINE
        I1=GRIDS(IMOD)%TRIGP(ITRI,1)
        I2=GRIDS(IMOD)%TRIGP(ITRI,2)
        I3=GRIDS(IMOD)%TRIGP(ITRI,3)
+
+       CALL FIX_PERIODCITY(I1,I2,I3,GRIDS(IMOD)%XYB,PT)
 ! coordinates of the first vertex A
-       x1=GRIDS(IMOD)%XYB(I1,1)
-       y1=GRIDS(IMOD)%XYB(I1,2)
+       x1 = PT(1,1)
+       y1 = PT(1,2)
 ! coordinates of the 2nd vertex B
-       x2=GRIDS(IMOD)%XYB(I2,1)
-       y2=GRIDS(IMOD)%XYB(I2,2)
+       x2 = PT(2,1)
+       y2 = PT(2,2)
 !coordinates of the 3rd vertex C
-       x3=GRIDS(IMOD)%XYB(I3,1)
-       y3=GRIDS(IMOD)%XYB(I3,2)
+       x3 = PT(3,1)
+       y3 = PT(3,2)
 !with M = (XTIN,YTIN) the target point ... 
 !vector product of AB and AC
        sg3=(y3-y1)*(x2-x1)-(x3-x1)*(y2-y1)
@@ -2729,4 +2739,95 @@ END SUBROUTINE
       RETURN    
       END SUBROUTINE SETUGIOBP
 !/ ------------------------------------------------------------------- /
+
+      SUBROUTINE FIX_PERIODCITY(I1,I2,I3,XYB,PT)
+!/
+!/                  +-----------------------------------+
+!/                  | WAVEWATCH III           NOAA/NCEP |
+!/                  |        Steven Brus                |
+!/                  |        Ali Abdolali               |
+!/                  |                        FORTRAN 90 |
+!/                  | Last update :         21-May-2020 |
+!/                  +-----------------------------------+
+!/
+!/    21-May-2020 : Origination.                        ( version 6.07 )
+!/    
+!/
+!  1. Purpose :
+!
+!     Adjust element longitude coordinates for elements straddling the
+!     dateline with distance of ~360 degrees
+!
+!  2. Method :
+!
+!     Detect if element has nodes on both sides of dateline and adjust
+!     coordinates so that all nodes have the same sign
+!
+!  3. Parameters :
+!
+!     Parameter list
+!     ----------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER, INTENT(IN) :: I1, I2, I3
+      DOUBLE PRECISION, INTENT(IN) :: XYB(:,:)
+      REAL*8, INTENT(OUT) :: PT(3,2)
+!     ----------------------------------------------------------------
+!
+!     Local variables.
+!     ----------------------------------------------------------------
+      INTEGER :: I
+      INTEGER :: R1GT180, R2GT180, R3GT180
+!     ----------------------------------------------------------------
+!
+!  4. Subroutines used :
+!
+
+!  5. Called by :
+!
+!      Name         Type  Module   Description
+!     ----------------------------------------------------------------
+!      SPATIAL_GRID Subr. W3TRIAM  Triangle area calculation
+!      NVECTRI      Subr. W3TRIAM  Edge length, angle, normal calcuation
+!      IS_IN_UNGRID Subr. W3TRIAM  Point in element calculation
+!     ----------------------------------------------------------------
+!
+!  6. Error messages :
+!
+!       None.
+!
+!  7. Remarks :
+!
+!  8. Structure :
+!
+!  9. Switches :
+!
+! 10. Source code :     
+!/ ------------------------------------------------------------------- /
+
+      PT(1,1) = XYB(I1,1)
+      PT(1,2) = XYB(I1,2)
+      PT(2,1) = XYB(I2,1)
+      PT(2,2) = XYB(I2,2)
+      PT(3,1) = XYB(I3,1)
+      PT(3,2) = XYB(I3,2)
+
+
+       R1GT180 = MERGE(1, 0, ABS(PT(3,1)-PT(2,1)).GT.180)
+       R2GT180 = MERGE(1, 0, ABS(PT(1,1)-PT(3,1)).GT.180)
+       R3GT180 = MERGE(1, 0, ABS(PT(2,1)-PT(1,1)).GT.180)
+ ! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the dateline
+ ! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole
+ ! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline
+
+
+      IF ( R1GT180 + R2GT180 == 2 ) THEN
+             PT(3,1)=PT(3,1)-SIGN(360.0,(PT(3,1)-PT(2,1)))
+      ELSE IF ( R2GT180 + R3GT180 == 2 ) THEN
+             PT(1,1)=PT(1,1)-SIGN(360.0,(PT(1,1)-PT(2,1)))
+      ELSE IF ( R1GT180 + R3GT180 == 2 ) THEN
+             PT(2,1)=PT(2,1)-SIGN(360.0,(PT(2,1)-PT(3,1)))
+      ENDIF
+
+      RETURN
+      END SUBROUTINE FIX_PERIODCITY
 END MODULE W3TRIAMD