Double precision option for regridded SMC output in ounf

model/inp/ww3_ounf.inp

@@ -26,10 +26,12 @@ $ netCDF version [3,4]
 $        and variable type 4 [2 = SHORT, 3 = it depends , 4 = REAL]
 $ swell partitions [0 1 2 3 4 5]
 $ variables in same file [T] or not [F]
+$ double precision output [T] or single precision [F] NOT CURRENTLY SUPPORTED IN .inp
 $
  3 4
  0 1 2
  T
+$ T
 $
 $ -------------------------------------------------------------------- $
 $ File prefix

model/nml/ww3_ounf.nml

@@ -99,6 +99,7 @@
 !     SMC%EYO           = -999.9         ! Last latitude
 !     SMC%CELFAC        = 1              ! Cell size factor (SMCTYPE=2 only)
 !     SMC%NOVAL         = UNDEF          ! Fill value for wet cells with no data
+!     SMC%DOUBLE_COORD               = F                  ! Return coordinates as double precision floats
 ! -------------------------------------------------------------------- !
 &SMC_NML
 /

model/src/w3gridmd.F90

@@ -632,10 +632,11 @@ MODULE W3GRIDMD
   !
   REAL                    :: RXFR, RFR1, SIGMA, SXFR, FACHF,      &
        VSC, VSC0, VOF,                      &
-       ZLIM, X, Y, XP,  XO0, YO0, DXO, DYO, &
+       ZLIM, X, Y, XP,                      &
        XO, YO, RD(4), RDTOT,                &
        FACTOR, RTH0, FMICHE, RWNDC,         &
        WCOR1, WCOR2
+  DOUBLE PRECISION        :: XO0, YO0, DXO, DYO
   !
   CHARACTER(LEN=4)        :: GSTRG, CSTRG
   !

model/src/w3metamd.F90

@@ -43,7 +43,7 @@ MODULE W3METAMD
   TYPE META_PAIR_T
     CHARACTER(LEN=64)  :: ATTNAME = UNSETC  !< Attribute name
     CHARACTER(LEN=120) :: ATTVAL = UNSETC   !< Attribute value
-    CHARACTER          :: TYPE = 'c'        !< Attribute type (c,i,f/r)
+    CHARACTER          :: TYPE = 'c'        !< Attribute type (c,i,f/r,d)
     TYPE(META_PAIR_T), POINTER  :: NEXT     !< Pointer to next node
   END TYPE META_PAIR_T
 
@@ -59,6 +59,7 @@ MODULE W3METAMD
     MODULE PROCEDURE META_LIST_APPEND_R !< Append a REAL value
     MODULE PROCEDURE META_LIST_APPEND_I !< Append an INTEGER value
     MODULE PROCEDURE META_LIST_APPEND_C !< Append a CHARACTER value
+    MODULE PROCEDURE META_LIST_APPEND_D !< Append a DOUBLE PRECISION value
   END INTERFACE META_LIST_APPEND
 
 CONTAINS
@@ -299,6 +300,35 @@ SUBROUTINE META_LIST_APPEND_C(LIST, ATTNAME, SVAL)
   END SUBROUTINE META_LIST_APPEND_C
 
 
+!/------------------------------------------------------------------ /
+!> @brief Append DOUBLE value attribute to list
+!>
+!> @param[in,out] LIST     The list to append to
+!> @param[in]     ATTNAME  The attribute name
+!> @param[in]     DVAL     The attribute value (DOUBLE)
+!>
+!> @author Kit Stokes
+!/------------------------------------------------------------------- /
+      SUBROUTINE META_LIST_APPEND_D(LIST, ATTNAME, DVAL)
+
+      IMPLICIT NONE
+
+      TYPE(META_LIST_T), INTENT(INOUT) :: LIST
+      CHARACTER(*), INTENT(IN)         :: ATTNAME
+      DOUBLE PRECISION, INTENT(IN)     :: DVAL
+!/------------------------------------------------------------------- /
+!/ Local parameters
+!/
+      TYPE(META_PAIR_T) :: META
+
+      META%ATTNAME = ATTNAME
+      WRITE(META%ATTVAL,*) DVAL
+      META%TYPE = 'd'
+      CALL META_LIST_APPEND(LIST, META)
+
+      END SUBROUTINE META_LIST_APPEND_D
+
+
   !/ ------------------------------------------------------------------- /
   !> @brief Find (first) entry in list with matching attname
   !>

model/src/w3nmlounfmd.F90

@@ -56,11 +56,12 @@ MODULE W3NMLOUNFMD
   ! smc grid structure
   TYPE NML_SMC_T
     INTEGER                     :: TYPE
-    REAL                        :: SXO
-    REAL                        :: SYO
+    DOUBLE PRECISION            :: SXO
+    DOUBLE PRECISION            :: SYO
     REAL                        :: EXO
     REAL                        :: EYO
     INTEGER                     :: CELFAC
+    LOGICAL                     :: DOUBLE_COORD
   END TYPE NML_SMC_T
 
   ! miscellaneous
@@ -515,6 +516,7 @@ SUBROUTINE READ_SMC_NML (NDSI, NML_SMC)
     SMC%EYO       = -999.9
     SMC%CELFAC    = 1
     SMC%TYPE      = 1
+    SMC%DOUBLE_COORD     = .FALSE.
 
     ! read smc namelist
     REWIND (NDSI)
@@ -806,8 +808,10 @@ SUBROUTINE REPORT_SMC_NML (NML_SMC)
     WRITE (NDSN,14) TRIM(MSG),'EXO       = ', NML_SMC%EXO
     WRITE (NDSN,14) TRIM(MSG),'EYO       = ', NML_SMC%EYO
     WRITE (NDSN,11) TRIM(MSG),'CELFAC    = ', NML_SMC%CELFAC
+    WRITE (NDSN,13) TRIM(MSG),'DOUBLE_COORD  = ', NML_SMC%DOUBLE_COORD
 
 11  FORMAT (A,2X,A,I12)
+13  FORMAT (A,2X,A,L1)
 14  FORMAT (A,2X,A,F8.2)
 
   END SUBROUTINE REPORT_SMC_NML

model/src/w3ounfmetamd.F90

@@ -1268,6 +1268,7 @@ SUBROUTINE GET_ATTVAL_TYPE(BUF, ILINE, ATTV, ATT_TYPE)
     !
     REAL :: R
     INTEGER :: I, IERR
+    DOUBLE PRECISION :: D
 
     ! Get attribute and type (default to "c" if no type set)
     ATT_TYPE = 'c'
@@ -1291,6 +1292,14 @@ SUBROUTINE GET_ATTVAL_TYPE(BUF, ILINE, ATTV, ATT_TYPE)
         WRITE(NDSE, 8001) "REAL/FLOAT", TRIM(FN_META), ILINE, TRIM(ATTV)
         CALL EXTCDE(10)
       ENDIF
+
+    CASE("d")
+      READ(attv, *, iostat=ierr) d
+      IF(ierr .ne. 0) THEN
+        WRITE(NDSE, 8001) "DOUBLE", TRIM(FN_META), ILINE, TRIM(ATTV)
+        CALL EXTCDE(10)
+      ENDIF
+
 
     CASE("c")
       ! Always ok.
@@ -1308,7 +1317,7 @@ SUBROUTINE GET_ATTVAL_TYPE(BUF, ILINE, ATTV, ATT_TYPE)
          '  => ', A /)
     !
 8002 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNFMETA : '/           &
-         '     ATTRIBUTE TYPE SHOULD BE ONE OF [c,i,r] '/       &
+         '     ATTRIBUTE TYPE SHOULD BE ONE OF [c,i,r,d] '/       &
          '     FILENAME = ', A /                                &
          '     LINE NO =', I5 /                                 &
          '  => ', A /)
@@ -1322,7 +1331,7 @@ END SUBROUTINE GET_ATTVAL_TYPE
   !>    or EOF is found. Splits meta pairs on the `=` character.
   !>
   !>    Freeform metadata pairs can also provide a variable type
-  !>    ("c", "i", or "r"; for character, int or real respectively).
+  !>    ("c", "i", "r", or "d"; for character, int, real, or double respectively).
   !>    String values with spaces should be quoted.
   !>
   !> @param[in]      NDMI      Unit number of metadata input file
@@ -2325,6 +2334,7 @@ SUBROUTINE WRITE_FREEFORM_META_LIST(NCID, VARID, METALIST, ERR)
     !/
     INTEGER :: I, IVAL
     REAL    :: RVAL
+    DOUBLE PRECISION :: DVAL
     TYPE(META_PAIR_T), POINTER :: P
 
     IF(METALIST%N .EQ. 0) RETURN
@@ -2343,6 +2353,11 @@ SUBROUTINE WRITE_FREEFORM_META_LIST(NCID, VARID, METALIST, ERR)
         READ(P%ATTVAL, *) IVAL
         ERR = NF90_PUT_ATT(NCID, VARID, P%ATTNAME, IVAL)
         IF(ERR /= NF90_NOERR) RETURN
+
+      CASE('d')
+        READ(P%ATTVAL, *) DVAL
+        ERR = NF90_PUT_ATT(NCID, VARID, P%ATTNAME, DVAL)
+        IF(ERR /= NF90_NOERR) RETURN
 
       CASE('r', 'f')
         READ(P%ATTVAL, *) RVAL

model/src/w3servmd.F90

@@ -64,6 +64,7 @@ MODULE W3SERVMD
   !     W3ACTURN       turns wave action(k,nth) anti-clockwise by AnglD.
   !     W3LLTOEQ       convert standard into rotated lat/lon, plus AnglD
   !     W3EQTOLL       revers of the LLTOEQ, but AnglD unchanged.
+  !     W3EQTOLL_DBL   double precision version of EQTOLL
   !     W3THTRN        turns direction value anti-clockwise by AnglD
   !     W3XYTRN        turns 2D vectors anti-clockwise by AnglD
   !
@@ -1658,6 +1659,116 @@ SUBROUTINE W3EQTOLL( PHI_EQ, LAMBDA_EQ, PHI, LAMBDA,   &
     RETURN
   END SUBROUTINE W3EQTOLL
 
+
+  ! Same as W3EQTOLL but working in double precision (D. Brown).
+  SUBROUTINE W3EQTOLL_DBL( PHI_EQ, LAMBDA_EQ, PHI, LAMBDA,   &
+       &                 ANGLED, PHI_POLE, LAMBDA_POLE, POINTS )
+
+    INTEGER:: POINTS      !IN  Number of points to be processed
+
+    DOUBLE PRECISION :: PHI_POLE,   & !IN  Latitude of equatorial lat-lon pole
+         &        LAMBDA_POLE   !IN  Longitude of equatorial lat-lon pole
+
+    DOUBLE PRECISION, DIMENSION(POINTS) ::         &
+         &        PHI,       & !OUT Latitude
+         &        LAMBDA,    & !OUT Longitude (0 =< LON < 360)
+         &        ANGLED,    & !OUT turning angle in deg for standard wind
+         &        LAMBDA_EQ, & !IN  Longitude in equatorial lat-lon coords
+         &        PHI_EQ       !IN  Latitude in equatorial lat-lon coords
+
+    ! Local varables:------------------------------------------------------
+    REAL(KIND=8) :: E_LAMBDA, E_PHI, A_LAMBDA, A_PHI,                 &
+         SIN_PHI_POLE, COS_PHI_POLE,                       &
+         TERM1, TERM2, ARG, LAMBDA_ZERO
+    INTEGER :: I
+
+    REAL(KIND=8), PARAMETER :: SMALL=1.0E-6
+
+    ! Double precision versions of values in constants.ftn:
+    REAL(KIND=8), PARAMETER         :: PI = 3.141592653589793
+    REAL(KIND=8), PARAMETER         :: RECIP_PI_OVER_180 = 180. / PI
+    REAL(KIND=8), PARAMETER         :: PI_OVER_180   = PI / 180.
+
+    ! ----------------------------------------------------------------------
+
+    ! 1. Initialise local constants
+
+    ! Latitude of zeroth meridian
+    LAMBDA_ZERO=LAMBDA_POLE+180.D0
+    ! Sine and cosine of latitude of eq pole
+    IF (PHI_POLE >= 0.0) THEN
+      SIN_PHI_POLE =  SIN(PI_OVER_180*PHI_POLE)
+      COS_PHI_POLE =  COS(PI_OVER_180*PHI_POLE)
+    ELSE
+      SIN_PHI_POLE = -SIN(PI_OVER_180*PHI_POLE)
+      COS_PHI_POLE = -COS(PI_OVER_180*PHI_POLE)
+    ENDIF
+
+    ! 2. Transform from equatorial to standard latitude-longitude
+
+    DO I= 1, POINTS
+
+      ! Scale eq longitude to range -180 to +180 degs
+
+      E_LAMBDA=LAMBDA_EQ(I)
+      IF(E_LAMBDA.GT. 180.0) E_LAMBDA=E_LAMBDA-360.D0
+      IF(E_LAMBDA.LT.-180.0) E_LAMBDA=E_LAMBDA+360.D0
+
+      ! Convert eq latitude & longitude to radians
+
+      E_LAMBDA=PI_OVER_180*E_LAMBDA
+      E_PHI=PI_OVER_180*PHI_EQ(I)
+
+      ! Compute latitude using equation (4.7)
+
+      ARG=COS_PHI_POLE*COS(E_PHI)*COS(E_LAMBDA) + SIN_PHI_POLE*SIN(E_PHI)
+      ARG=MIN(ARG, 1.D0)
+      ARG=MAX(ARG,-1.D0)
+      A_PHI=ASIN(ARG)
+      PHI(I)=RECIP_PI_OVER_180*A_PHI
+
+      ! Compute longitude using equation (4.8)
+
+      TERM1 = COS(E_PHI)*SIN_PHI_POLE*COS(E_LAMBDA) - SIN(E_PHI)*COS_PHI_POLE
+      TERM2 = COS(A_PHI)
+      IF(TERM2.LT.SMALL) THEN
+        A_LAMBDA=0.D0
+      ELSE
+        ARG=TERM1/TERM2
+        ARG=MIN(ARG, 1.D0)
+        ARG=MAX(ARG,-1.D0)
+        A_LAMBDA=RECIP_PI_OVER_180*ACOS(ARG)
+        A_LAMBDA=SIGN(A_LAMBDA,E_LAMBDA)
+        A_LAMBDA=A_LAMBDA+LAMBDA_ZERO
+      END IF
+
+      ! Scale longitude to range 0 to 360 degs
+
+      IF(A_LAMBDA.GE.360.0) A_LAMBDA=A_LAMBDA-360.D0
+      IF(A_LAMBDA.LT.  0.0) A_LAMBDA=A_LAMBDA+360.D0
+      LAMBDA(I)=A_LAMBDA
+
+      !Li  Calculate turning angle for standard wind velocity
+
+      A_LAMBDA=PI_OVER_180*(LAMBDA(I)-LAMBDA_ZERO)
+
+      ! Formulae used are from eqs (4.19) and (4.21)
+
+      TERM2=SIN(E_LAMBDA)
+      ARG=SIN(A_LAMBDA)*TERM2*SIN_PHI_POLE     &
+           &           +COS(A_LAMBDA)*COS(E_LAMBDA)
+      ARG=MIN(ARG, 1.D0)
+      ARG=MAX(ARG,-1.D0)
+      TERM1=RECIP_PI_OVER_180*ACOS(ARG)
+      ANGLED(I)=SIGN(TERM1,TERM2)
+      !Li
+
+    ENDDO
+
+    RETURN
+  END SUBROUTINE W3EQTOLL_DBL
+
+
   !Li
   !/ ------------------------------------------------------------------- /
   !/ ------------------------------------------------------------------- /

model/src/w3smcomd.F90

@@ -587,7 +587,7 @@ SUBROUTINE READ_SMCINT()
     ! output to a rotated pole grid...
 
     NDSMC = 50
-    OPEN(NDSMC, file='smcint.ww3', status='old', form='unformatted', convert=file_endian, iostat=ierr)
+    OPEN(NDSMC, file='smcint.ww3', status='old', form='unformatted', convert='native', iostat=ierr)
     IF(ierr .NE. 0) THEN
       WRITE(*,*) "ERROR! Failed to open smcint.ww3 for reading"
       CALL EXTCDE(1)