From 8b5e91f0887ab2b8c218228fff6cb82297a613b3 Mon Sep 17 00:00:00 2001 From: Mickael Accensi <49198861+mickaelaccensi@users.noreply.github.com> Date: Thu, 4 Apr 2024 21:34:11 +0200 Subject: [PATCH] add output parameters for skewness (#1209) Co-authored-by: Fabrice Ardhuin --- manual/eqs/output.tex | 4 + manual/manual.bib | 10 + model/inp/ww3_ounf.inp | 2 +- model/inp/ww3_shel.inp | 11 +- model/nml/ww3_ounf.nml | 2 +- model/nml/ww3_shel.nml | 3 + model/src/w3adatmd.F90 | 42 +- model/src/w3initmd.F90 | 86 ++++- model/src/w3iogomd.F90 | 556 ++++++++++++++++++++++++++- model/src/w3odatmd.F90 | 5 +- model/src/w3ounfmetamd.F90 | 39 ++ model/src/ww3_ounf.F90 | 15 +- model/src/ww3_outf.F90 | 34 ++ model/tools/bash/ww3_ounf_inp2nml.sh | 2 +- model/tools/bash/ww3_shel_inp2nml.sh | 3 + regtests/ww3_ts1/input/ww3_ounf.inp | 2 +- regtests/ww3_ts1/input/ww3_ounf.nml | 2 +- regtests/ww3_ts1/input/ww3_shel.inp | 2 +- regtests/ww3_ts1/input/ww3_shel.nml | 2 +- 19 files changed, 800 insertions(+), 22 deletions(-) diff --git a/manual/eqs/output.tex b/manual/eqs/output.tex index bfa7e0b5a..8deb484e1 100644 --- a/manual/eqs/output.tex +++ b/manual/eqs/output.tex @@ -310,6 +310,10 @@ \subsection{~Output parameters} \label{sub:outpars} \begin{equation} Q_{kk} = \frac{1}{E^2} \int_0^{f_{NK}} \int_0^{2\pi} 0.5 \left[ A(k,\theta)+ A(k,\theta+\pi)\right]^2 \frac{\sigma^2}{k C_g} \:\rd \theta \: \rd \sigma \: \label{eq:qkk} \end{equation} +\item \textbf{SKW} Skewness of surface elevation sampled at zero slope. This is the $\lambda_1$ parameter defined in \cite{Barrick&Lipa1985} or $\lambda_{3,0,0}$ in \cite{Srokosz1986}. It is computed from the second order correction to the surface elevation, using ECWAM code by P. Janssen. +\item \textbf{EMB} this is $-\gamma/8 = -(\lambda_{1,2,0}+\lambda_{1,0,2}-2 \lambda{0,1,1} \lambda{1,1,1})/8 (1-\lambda_{0,1,1]^2)$, such that the mean sea level of points with zero slope +is EMB$\times H_s$. +\item \textbf{EMC} this is hte additional tracker bias coefficient equal to $-\lambda_{3,0,0}/24$, which is specific to the choice of retracker, see the $J_z$ function in \cite{DeCarlo&Ardhuin2024}. \end{list} \item{Numerical diagnostics } diff --git a/manual/manual.bib b/manual/manual.bib index 3da650ea2..33e0a9fdd 100644 --- a/manual/manual.bib +++ b/manual/manual.bib @@ -3760,3 +3760,13 @@ @PHDTHESIS{Gagnaire-Renou2009 year = 2010, } +@ARTICLE{Srokosz1986, + author = "Meric A. Srokosz", + title = "On the joint distribution of surface elevation and slopes for a non linear random sea, with an application to radar altimetry", + journal = JGR, + volume = 91, + pages = "995--1006", + year = 1986, + keywords={altimeter;sea state bias}, +} + diff --git a/model/inp/ww3_ounf.inp b/model/inp/ww3_ounf.inp index 7bde30754..b0c29ff75 100644 --- a/model/inp/ww3_ounf.inp +++ b/model/inp/ww3_ounf.inp @@ -16,7 +16,7 @@ $ DPT CUR WND AST WLV ICE IBG TAU RHO D50 IC1 IC5 HS LM T02 T0M1 T01 FP $ DIR SPR DP HIG EF TH1M STH1M TH2M STH2M WN PHS PTP PLP PDIR PSPR PWS PDP $ PQP PPE PGW PSW PTM10 PT01 PT02 PEP TWS PNR UST CHA CGE FAW TAW TWA WCC $ WCF WCH WCM SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC ABR UBR BED -$ FBB TBB MSS MSC DTD FC CFX CFD CFK U1 U2 WNM TOC +$ FBB TBB MSS MSC DTD FC CFX CFD CFK U1 U2 WNM TOC MSS QP QKK SKW EMB EMC $ N DPT HS FP T01 diff --git a/model/inp/ww3_shel.inp b/model/inp/ww3_shel.inp index 7980cbe06..576e01c5f 100644 --- a/model/inp/ww3_shel.inp +++ b/model/inp/ww3_shel.inp @@ -212,10 +212,13 @@ $ 8 Spectrum parameters $ ------------------------------------------------- $ F F 8 1 MSS[X,Y] MSS Mean square slopes $ F F 8 2 MSC[X,Y] MSC Spectral level at high frequency tail -! F F 8 3 MSSD MSD Slope direction -! F F 8 4 MSCD MCD Tail slope direction -! F F 8 5 QP QP Goda peakedness parameter -! F F 8 6 QKK QKK Wavenumber peakedness +$ F F 8 3 MSSD MSD Slope direction +$ F F 8 4 MSCD MCD Tail slope direction +$ F F 8 5 QP QP Goda peakedness parameter +$ F F 8 6 QKK QKK Wavenumber peakedness +$ F F 8 7 SKEW SKW Skewness of elevation for zero slopes +$ F F 8 8 EMBIA1 EMB Mean sea level at zero slopes / Hs +$ F F 8 9 EMBIA2 EMC Tracker bias for LRM least square altimetry $ ------------------------------------------------- $ 9 Numerical diagnostics $ ------------------------------------------------- diff --git a/model/nml/ww3_ounf.nml b/model/nml/ww3_ounf.nml index 9b1ffe136..a70cf9eae 100644 --- a/model/nml/ww3_ounf.nml +++ b/model/nml/ww3_ounf.nml @@ -14,7 +14,7 @@ ! UST CHA CGE FAW TAW TWA WCC WCF WCH WCM FWS ! SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC TOC ! ABR UBR BED FBB TBB -! MSS MSC WL02 AXT AYT AXY +! MSS MSC MSD MCD QP QKK SKW EMB EMC ! DTD FC CFX CFD CFK ! U1 U2 ! diff --git a/model/nml/ww3_shel.nml b/model/nml/ww3_shel.nml index 97beaf6a0..b6ae8cfef 100644 --- a/model/nml/ww3_shel.nml +++ b/model/nml/ww3_shel.nml @@ -206,6 +206,9 @@ ! F F 8 4 MSCD MCD Tail slope direction ! F F 8 5 QP QP Goda peakedness parameter ! F F 8 6 QKK QKK Wavenumber peakedness +! F F 8 7 SKEW SKW Skewness of elevation for zero slopes +! F F 8 8 EMBIA1 EMB Mean sea level at zero slopes / Hs +! F F 8 9 EMBIA2 EMC Tracker bias for LRM least square altimetry ! ------------------------------------------------- ! 9 Numerical diagnostics ! ------------------------------------------------- diff --git a/model/src/w3adatmd.F90 b/model/src/w3adatmd.F90 index 1ee07eca4..2daee3609 100644 --- a/model/src/w3adatmd.F90 +++ b/model/src/w3adatmd.F90 @@ -188,6 +188,7 @@ MODULE W3ADATMD ! MSCD R.A. Public Direction of MSCX ! QP R.A. Public Goda peakedness parameter. ! QKK R.A. Public Spectral bandwidth (De Carlo et al. 2023) + ! SKEW R.A. Public skewness lambda_3,0,0 (Srokosz 1986) ! ! DTDYN R.A. Public Mean dynamic time step (raw). ! FCUT R.A. Public Cut-off frequency for tail. @@ -475,9 +476,10 @@ MODULE W3ADATMD ! Output fields group 8) ! REAL, POINTER :: MSSX(:), MSSY(:), MSSD(:), & - MSCX(:), MSCY(:), MSCD(:), QKK(:) + MSCX(:), MSCY(:), MSCD(:), QKK(:), SKEW(:), EMBIA1(:), EMBIA2(:) REAL, POINTER :: XMSSX(:), XMSSY(:), XMSSD(:), & - XMSCX(:), XMSCY(:), XMSCD(:), XQKK(:) + XMSCX(:), XMSCY(:), XMSCD(:), XQKK(:), & + XSKEW(:), XEMBIA1(:), XEMBIA2(:) ! ! Output fields group 9) ! @@ -613,7 +615,7 @@ MODULE W3ADATMD BEDFORMS(:,:), PHIBBL(:), TAUBBL(:,:) ! REAL, POINTER :: MSSX(:), MSSY(:), MSSD(:), & - MSCX(:), MSCY(:), MSCD(:), QKK(:) + MSCX(:), MSCY(:), MSCD(:), QKK(:), SKEW(:), EMBIA1(:), EMBIA2(:) ! REAL, POINTER :: DTDYN(:), FCUT(:), CFLXYMAX(:), & CFLTHMAX(:), CFLKMAX(:) @@ -1265,7 +1267,9 @@ SUBROUTINE W3DIMA ( IMOD, NDSE, NDST, D_ONLY ) ALLOCATE ( WADATS(IMOD)%MSSX(NSEALM), WADATS(IMOD)%MSSY(NSEALM), & WADATS(IMOD)%MSCX(NSEALM), WADATS(IMOD)%MSCY(NSEALM), & WADATS(IMOD)%MSSD(NSEALM), WADATS(IMOD)%MSCD(NSEALM), & - WADATS(IMOD)%QKK(NSEALM), STAT=ISTAT ) + WADATS(IMOD)%QKK(NSEALM), WADATS(IMOD)%SKEW(NSEALM), & + WADATS(IMOD)%EMBIA1(NSEALM), WADATS(IMOD)%EMBIA2(NSEALM), & + STAT=ISTAT ) CHECK_ALLOC_STATUS ( ISTAT ) ! WADATS(IMOD)%MSSX = UNDEF @@ -1275,6 +1279,9 @@ SUBROUTINE W3DIMA ( IMOD, NDSE, NDST, D_ONLY ) WADATS(IMOD)%MSCY = UNDEF WADATS(IMOD)%MSCD = UNDEF WADATS(IMOD)%QKK = UNDEF + WADATS(IMOD)%SKEW = UNDEF + WADATS(IMOD)%EMBIA1 = UNDEF + WADATS(IMOD)%EMBIA2 = UNDEF call print_memcheck(memunit, 'memcheck_____:'//' W3DIMA 8') ! ! 9) Numerical diagnostics @@ -2281,6 +2288,24 @@ SUBROUTINE W3XDMA ( IMOD, NDSE, NDST, OUTFLAGS ) ALLOCATE ( WADATS(IMOD)%XQKK(1) ) END IF ! + IF ( OUTFLAGS( 8, 7) ) THEN + ALLOCATE ( WADATS(IMOD)%XSKEW(NXXX) ) + ELSE + ALLOCATE ( WADATS(IMOD)%XSKEW(1) ) + END IF + ! + IF ( OUTFLAGS( 8, 8) ) THEN + ALLOCATE ( WADATS(IMOD)%XEMBIA1(NXXX) ) + ELSE + ALLOCATE ( WADATS(IMOD)%XEMBIA1(1) ) + END IF + ! + IF ( OUTFLAGS( 8, 9) ) THEN + ALLOCATE ( WADATS(IMOD)%XEMBIA2(NXXX) ) + ELSE + ALLOCATE ( WADATS(IMOD)%XEMBIA2(1) ) + END IF + ! WADATS(IMOD)%XMSSX = UNDEF WADATS(IMOD)%XMSSY = UNDEF WADATS(IMOD)%XMSSD = UNDEF @@ -2289,6 +2314,9 @@ SUBROUTINE W3XDMA ( IMOD, NDSE, NDST, OUTFLAGS ) WADATS(IMOD)%XMSCD = UNDEF WADATS(IMOD)%XQP(1) = UNDEF WADATS(IMOD)%XQKK = UNDEF + WADATS(IMOD)%XSKEW = UNDEF + WADATS(IMOD)%XEMBIA1 = UNDEF + WADATS(IMOD)%XEMBIA2 = UNDEF ! IF ( OUTFLAGS( 9, 1) ) THEN ALLOCATE ( WADATS(IMOD)%XDTDYN(NXXX), STAT=ISTAT ) @@ -2903,6 +2931,9 @@ SUBROUTINE W3SETA ( IMOD, NDSE, NDST ) MSCY => WADATS(IMOD)%MSCY MSCD => WADATS(IMOD)%MSCD QKK => WADATS(IMOD)%QKK + SKEW => WADATS(IMOD)%SKEW + EMBIA1 => WADATS(IMOD)%EMBIA1 + EMBIA2 => WADATS(IMOD)%EMBIA2 ! DTDYN => WADATS(IMOD)%DTDYN FCUT => WADATS(IMOD)%FCUT @@ -3242,6 +3273,9 @@ SUBROUTINE W3XETA ( IMOD, NDSE, NDST ) MSCY => WADATS(IMOD)%XMSCY MSCD => WADATS(IMOD)%XMSCD QKK => WADATS(IMOD)%XQKK + SKEW => WADATS(IMOD)%XSKEW + EMBIA1 => WADATS(IMOD)%XEMBIA1 + EMBIA2 => WADATS(IMOD)%XEMBIA2 ! DTDYN => WADATS(IMOD)%XDTDYN FCUT => WADATS(IMOD)%XFCUT diff --git a/model/src/w3initmd.F90 b/model/src/w3initmd.F90 index 93218d473..044a18760 100644 --- a/model/src/w3initmd.F90 +++ b/model/src/w3initmd.F90 @@ -2150,7 +2150,7 @@ SUBROUTINE W3MPIO ( IMOD ) STMAXE, STMAXD, HMAXE, HCMAXE, HMAXD, & HCMAXD, QP, PTHP0, PQP, PPE, PGW, PSW, & PTM1, PT1, PT2, PEP, WBT, CX, CY, & - TAUOCX, TAUOCY, WNMEAN, QKK + TAUOCX, TAUOCY, WNMEAN, QKK, SKEW, EMBIA1, EMBIA2 #endif #ifdef W3_MPI @@ -3406,6 +3406,48 @@ SUBROUTINE W3MPIO ( IMOD ) #ifdef W3_MPIT WRITE (NDST,9011) IH, ' 8/06', IROOT, IT, IRQGO(IH), IERR #endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 7) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_SEND_INIT (SKEW (1),NSEALM , MPI_REAL, IROOT, & + IT, MPI_COMM_WAVE, IRQGO(IH), IERR) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/07', IROOT, IT, IRQGO(IH), IERR +#endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 8) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_SEND_INIT (EMBIA1 (1),NSEALM , MPI_REAL, IROOT, & + IT, MPI_COMM_WAVE, IRQGO(IH), IERR) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/08', IROOT, IT, IRQGO(IH), IERR +#endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 9) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_SEND_INIT (EMBIA2 (1),NSEALM , MPI_REAL, IROOT, & + IT, MPI_COMM_WAVE, IRQGO(IH), IERR) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/09', IROOT, IT, IRQGO(IH), IERR +#endif #ifdef W3_MPI END IF #endif @@ -4653,6 +4695,48 @@ SUBROUTINE W3MPIO ( IMOD ) #ifdef W3_MPIT WRITE (NDST,9011) IH, ' 8/06', IFROM, IT, IRQGO2(IH), IERR #endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 7) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_RECV_INIT (SKEW (I0),1,WW3_FIELD_VEC, IFROM, IT, & + MPI_COMM_WAVE, IRQGO2(IH), IERR ) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/07', IFROM, IT, IRQGO2(IH), IERR +#endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 8) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_RECV_INIT (EMBIA1 (I0),1,WW3_FIELD_VEC, IFROM, IT, & + MPI_COMM_WAVE, IRQGO2(IH), IERR ) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/08', IFROM, IT, IRQGO2(IH), IERR +#endif +#ifdef W3_MPI + END IF +#endif + ! +#ifdef W3_MPI + IF ( FLGRDALL( 8, 9) ) THEN + IH = IH + 1 + IT = IT + 1 + CALL MPI_RECV_INIT (EMBIA2 (I0),1,WW3_FIELD_VEC, IFROM, IT, & + MPI_COMM_WAVE, IRQGO2(IH), IERR ) +#endif +#ifdef W3_MPIT + WRITE (NDST,9011) IH, ' 8/09', IFROM, IT, IRQGO2(IH), IERR +#endif #ifdef W3_MPI END IF #endif diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index 2ddfa77e0..de660ded4 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -18,7 +18,7 @@ MODULE W3IOGOMD !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | - !/ | Last update : 22-Mar-2021 | + !/ | Last update : 02-Mar-2024 | !/ +-----------------------------------+ !/ !/ 04-Jan-2001 : Origination. ( version 2.00 ) @@ -74,8 +74,9 @@ MODULE W3IOGOMD !/ 22-Mar-2021 : Add extra coupling fields as output ( version 7.13 ) !/ 21-Jul-2022 : Correct FP0 calc for peak energy in ( version 7.14 ) !/ min/max freq band (B. Pouliot, CMC) + !/ 02-Mar-2024 : Add skweness and EM bias varaible ( version 7.xx ) !/ - !/ Copyright 2009-2014 National Weather Service (NWS), + !/ Copyright 2009-2024 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. @@ -1126,6 +1127,15 @@ SUBROUTINE W3FLDTOIJ(FLD, I, J, IAPROC, NAPOUT, NDSEN) CASE('QKK') I = 8 J = 6 + CASE('SKW') + I = 8 + J = 7 + CASE('EMB') + I = 8 + J = 8 + CASE('EMC') + I = 8 + J = 9 ! ! Group 9 ! @@ -2340,6 +2350,11 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) IF (FLOLOC( 6, 12)) THEN CALL CALC_U3STOKES(A,2) ENDIF + ! + IF (FLOLOC( 8, 7).OR.FLOLOC( 8, 8).OR.FLOLOC( 8, 9)) THEN + CALL SKEWNESS(A) + END IF + ! ! Dominant wave breaking probability ! @@ -2417,6 +2432,7 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & !/ processing code) !/ 25-Aug-2018 : Add WBT parameter ( version 6.06 ) !/ 22-Mar-2021 : Add extra coupling fields as output ( version 7.13 ) + !/ 07-Mar-2024 : Add Skewness parameters ( version 7.13 ) !/ ! 1. Purpose : ! @@ -2513,7 +2529,7 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & CFLXYMAX, CFLTHMAX, CFLKMAX, P2SMS, US3D, & TH1M, STH1M, TH2M, STH2M, HSIG, PHICE, TAUICE,& STMAXE, STMAXD, HMAXE, HCMAXE, HMAXD, HCMAXD,& - USSP, TAUOCX, TAUOCY, QKK + USSP, TAUOCX, TAUOCY, QKK, SKEW, EMBIA1, EMBIA2 !/ USE W3ODATMD, ONLY: NOGRP, NGRPP, IDOUT, UNDEF, NDST, NDSE, & FLOGRD, IPASS => IPASS1, WRITE => WRITE1, & @@ -2923,6 +2939,9 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & IF ( FLOGRD( 8, 4) ) MSCD (ISEA) = UNDEF IF ( FLOGRD( 8, 5) ) QP (ISEA) = UNDEF IF ( FLOGRD( 8, 6) ) QKK (ISEA) = UNDEF + IF ( FLOGRD( 8, 7) ) SKEW (ISEA) = UNDEF + IF ( FLOGRD( 8, 8) ) EMBIA1(ISEA) = UNDEF + IF ( FLOGRD( 8, 9) ) EMBIA2(ISEA) = UNDEF ! IF ( FLOGRD( 9, 1) ) DTDYN (ISEA) = UNDEF IF ( FLOGRD( 9, 2) ) FCUT (ISEA) = UNDEF @@ -3576,7 +3595,7 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & #endif ! ! Section 8) - ! + !Skewness ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 1 ) THEN WRITE ( NDSOG ) MSSX(1:NSEA) #ifdef W3_ASCII @@ -3614,6 +3633,21 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & WRITE ( NDSOG ) QKK(1:NSEA) #ifdef W3_ASCII WRITE ( NDSOA,* ) 'QKK:', QKK(1:NSEA) +#endif + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 7 ) THEN + WRITE ( NDSOG ) SKEW(1:NSEA) +#ifdef W3_ASCII + WRITE ( NDSOA,* ) 'SKW:', SKEW(1:NSEA) +#endif + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 8 ) THEN + WRITE ( NDSOG ) EMBIA1(1:NSEA) +#ifdef W3_ASCII + WRITE ( NDSOA,* ) 'EMB:', EMBIA1(1:NSEA) +#endif + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 9 ) THEN + WRITE ( NDSOG ) EMBIA2(1:NSEA) +#ifdef W3_ASCII + WRITE ( NDSOA,* ) 'EMC:', EMBIA2(1:NSEA) #endif ! ! Section 9) @@ -3967,6 +4001,12 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD & READ (NDSOG,END=801,ERR=802,IOSTAT=IERR) QP(1:NSEA) ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 6 ) THEN READ (NDSOG,END=801,ERR=802,IOSTAT=IERR) QKK(1:NSEA) + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 7 ) THEN + READ (NDSOG,END=801,ERR=802,IOSTAT=IERR) SKEW(1:NSEA) + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 8 ) THEN + READ (NDSOG,END=801,ERR=802,IOSTAT=IERR) EMBIA1(1:NSEA) + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 9 ) THEN + READ (NDSOG,END=801,ERR=802,IOSTAT=IERR) EMBIA2(1:NSEA) ! ! Section 9) ! @@ -4610,4 +4650,512 @@ SUBROUTINE CALC_WBT (A) !/ END SUBROUTINE CALC_WBT !/ ------------------------------------------------------------------- / + !/ + !> + !> @brief Computation of second order harmonics and + !> relevant tables for the altimeter corrections + !> + !> @param[in] NKHF Extended number of frequencies. + !> @param[out] FAC0 2nd order coef correction. + !> @param[out] FAC1 2nd order coef correction. + !> @param[out] FAC2 2nd order coef correction. + !> @param[out] FAC3 2nd order coef correction. + !> + !> @author P. Janssen @date 29-Mar-2024 + !> + SUBROUTINE SECONDHH(NKHF,FAC0,FAC1,FAC2,FAC3) +!---------------------------------------------------------------- + +!**** *SECONDHH* - COMPUTATION OF SECOND ORDER HARMONICS AND +! RELEVANT TABLES FOR THE ALTIMETER CORRECTIONS. + +! P.A.E.M. JANSSEN + +! PURPOSE. +! --------- + +! COMPUTE SECOND HARMONICS + +!** INTERFACE. +! ---------- + +! *CALL* *SECONDHH* + +! METHOD. +! ------- + +! SEE REFERENCE. + +! EXTERNALS. +! ---------- + +! VMIN_D +! VPLUS_D + +! REFERENCES. +! ----------- + +! V E ZAKHAROV(1967) + +!------------------------------------------------------------------- + +!------------------------------------------------------------------- +USE CONSTANTS, ONLY: GRAV, TPI +USE W3GDATMD, ONLY: NK, NTH, XFR, SIG, TH, DTH, ECOS, ESIN + IMPLICIT NONE + ! REAL(KIND=4) :: VMIN_D,VPLUS_D + + + + INTEGER, INTENT(IN) :: NKHF + REAL(KIND=4), DIMENSION(NTH,NTH,NKHF,NKHF), INTENT(OUT) :: FAC0, FAC1, FAC2, FAC3 + REAL(KIND=4), PARAMETER :: FRATIO = 1.1 + + + INTEGER :: M, K1, M1, K2, M2 + + REAL(KIND=4), PARAMETER :: DEL1=1.0E-8 + REAL(KIND=4), PARAMETER :: ZCONST = 0.0281349 + + !REAL(KIND=4) :: VMIN_D, VPLUS_D + REAL(KIND=4) :: CO1 + REAL(KIND=4) :: XK1, XK1SQ, XK2, XK2SQ, XK3 + REAL(KIND=4) :: COSDIFF + REAL(KIND=4) :: X12, X13, X32, OM1, OM2, OM3, F1, F2, F3 + REAL(KIND=4) :: VM, VP + REAL(KIND=4) :: DELOM1, DELOM2 + REAL(KIND=4) :: DELOM321, DELOM312 + REAL(KIND=4) :: C22, S22 + + REAL(KIND=4), DIMENSION(NTH,NTH,NKHF,NKHF) :: B + REAL(KIND=4), DIMENSION(:), ALLOCATABLE:: FAK, SIGHF, DFIMHF + + + + +!----------------------------------------------------------------------- + + + + +!* 1. INITIALISE RELEVANT QUANTITIES. + + ALLOCATE(FAK(NKHF)) + ALLOCATE(SIGHF(NKHF)) + ALLOCATE(DFIMHF(NKHF)) + + SIGHF(1) = SIG(1) + DO M=2,NKHF + SIGHF(M) = XFR*SIGHF(M-1) + ENDDO + + DO M=1,NKHF + FAK(M) = (SIGHF(M))**2/GRAV + ENDDO + + CO1 = 0.5*(XFR-1.)*DTH + DFIMHF(1) = CO1*SIGHF(1) + DO M=2,NKHF-1 + DFIMHF(M)=CO1*(SIGHF(M)+SIGHF(M-1)) + ENDDO + DFIMHF(NKHF)=CO1*SIGHF(NKHF-1) + + DO M2=1,NKHF + XK2 = FAK(M2) + XK2SQ = FAK(M2)**2 + DO M1=1,NKHF + XK1 = FAK(M1) + XK1SQ = FAK(M1)**2 + DO K1=1,NTH + DO K2=1,NTH + COSDIFF = COS(TH(K1)-TH(K2)) + X12 = XK1*XK2*COSDIFF + XK3 = XK1SQ + XK2SQ +2.0*X12 +DEL1 + XK3 = SQRT(XK3) + X13 = XK1SQ+X12 + X32 = X12+XK2SQ + OM1 = SQRT(GRAV*XK1) + OM2 = SQRT(GRAV*XK2) + OM3 = SQRT(GRAV*XK3) + F1 = SQRT(XK1/(2.0*OM1)) + F2 = SQRT(XK2/(2.0*OM2)) + F3 = SQRT(XK3/(2.0*OM3)) + VM = TPI*VMIN_D(XK3,XK1,XK2,X13,X32,X12,OM3,OM1,OM2) + VP = TPI*VPLUS_D(-XK3,XK1,XK2,-X13,-X32,X12,OM3,OM1,OM2) + DELOM1 = OM3-OM1-OM2+DEL1 + DELOM2 = OM3+OM1+OM2+DEL1 + FAC0(K1,K2,M1,M2) = -F3/(F1*F2)*(VM/(DELOM1)+ & + & VP/(DELOM2)) + ENDDO + ENDDO + ENDDO + ENDDO + + DO M2=1,NKHF + XK2 = FAK(M2) + XK2SQ = FAK(M2)**2 + DO M1=1,NKHF + XK1 = FAK(M1) + XK1SQ = FAK(M1)**2 + DO K1=1,NTH + DO K2=1,NTH + COSDIFF = COS(TH(K1)-TH(K2)) + X12 = XK1*XK2*COSDIFF + XK3 = XK1SQ + XK2SQ - 2.*X12 + DEL1 + XK3 = SQRT(XK3) + X13 = XK1SQ-X12 + X32 = X12-XK2SQ + OM1 = SQRT(GRAV*XK1) + OM2 = SQRT(GRAV*XK2) + OM3 = SQRT(GRAV*XK3)+DEL1 + F1 = SQRT(XK1/(2.0*OM1)) + F2 = SQRT(XK2/(2.0*OM2)) + F3 = SQRT(ABS(XK3)/(2.0*OM3)) + VM = TPI*VMIN_D(XK1,XK3,XK2,X13,X12,X32,OM1,OM3,OM2) + VP = TPI*VMIN_D(XK2,-XK3,XK1,-X32,X12,-X13,OM2,OM3,OM1) + DELOM321 = OM3+OM2-OM1+DEL1 + DELOM312 = OM3+OM1-OM2+DEL1 + B(K1,K2,M1,M2) = -F3/(F1*F2)*(VM/(DELOM321)+ & + & VP/(DELOM312)) + ENDDO + ENDDO + ENDDO + ENDDO + + DO M2=1,NKHF + XK2SQ = FAK(M2)**2 + DO M1=1,NKHF + XK1SQ = FAK(M1)**2 + DO K2=1,NTH + DO K1=1,NTH + C22 = FAC0(K1,K2,M1,M2)+B(K1,K2,M1,M2) + S22 = B(K1,K2,M1,M2)-FAC0(K1,K2,M1,M2) + FAC1(K1,K2,M1,M2) = & + & (XK1SQ*ECOS(K1)**2 + XK2SQ*ECOS(K2)**2)*C22 & + & -FAK(M1)*FAK(M2)*ECOS(K1)*ECOS(K2)*S22 + FAC2(K1,K2,M1,M2) = & + & (XK1SQ*ESIN(K1)**2 + XK2SQ*ESIN(K2)**2)*C22 & + & -FAK(M1)*FAK(M2)*ESIN(K1)*ESIN(K2)*S22 + FAC3(K1,K2,M1,M2) = & + & (XK1SQ*ESIN(K1)*ECOS(K1) + & + & XK2SQ*ESIN(K2)*ECOS(K2))*C22 & + & -FAK(M1)*FAK(M2)*ECOS(K1)*ESIN(K2)*S22 + FAC0(K1,K2,M1,M2) = C22 + ENDDO + ENDDO + ENDDO + ENDDO + + + CONTAINS + +!----------------------------------------------------------------------- + + REAL(KIND=4) FUNCTION VMIN_D(XI,XJ,XK,XIJ,XIK,XJK,XOI,XOJ,XOK) + +! PETER JANSSEN + +! PURPOSE. +! -------- + +! GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE +! WAVE INTERACTIONS OF DEEP-WATER WAVES IN THE +! IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV) + +! INTERFACE. +! ---------- +! *VMIN_D(XI,XJ,XK)* +! *XI* - WAVE NUMBER +! *XJ* - WAVE NUMBER +! *XK* - WAVE NUMBER +! METHOD. +! ------- +! NONE + +! EXTERNALS. +! ---------- +! NONE. + + +!*** 1. DETERMINE NONLINEAR TRANSFER. +! -------------------------------- + IMPLICIT NONE + REAL, INTENT(IN) :: XI, XJ, XK, XIJ, XIK, XJK, XOI, XOJ, XOK + REAL :: RI, RJ, RK, OI, OJ, OK, SQIJK, SQIKJ, SQJKI + + RI=ABS(XI)+DEL1 + RJ=ABS(XJ)+DEL1 + RK=ABS(XK)+DEL1 + OI=XOI+DEL1 + OJ=XOJ+DEL1 + OK=XOK+DEL1 + SQIJK=SQRT(OI*OJ*RK/(OK*RI*RJ)) + SQIKJ=SQRT(OI*OK*RJ/(OJ*RI*RK)) + SQJKI=SQRT(OJ*OK*RI/(OI*RJ*RK)) + VMIN_D=ZCONST*( (XIJ-RI*RJ)*SQIJK + (XIK-RI*RK)*SQIKJ & + & + (XJK+RJ*RK)*SQJKI ) + + END FUNCTION VMIN_D + +!----------------------------------------------------------------------- + + REAL(KIND=4) FUNCTION VPLUS_D(XI,XJ,XK,XIJ,XIK,XJK,XOI,XOJ,XOK) + +!*** *VPLUS_D* DETERMINES THE NONLINEAR TRANSFER COEFFICIENT FOR THREE +! WAVE INTERACTIONS OF DEEP-WATER WAVES. + +! PETER JANSSEN + +! PURPOSE. +! -------- + +! GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE +! WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE +! IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV) + +! INTERFACE. +! ---------- +! *VPLUS_D(XI,XJ,XK)* +! *XI* - WAVE NUMBER +! *XJ* - WAVE NUMBER +! *XK* - WAVE NUMBER +! METHOD. +! ------- +! NONE + +! EXTERNALS. +! ---------- +! NONE. + + + +!*** 1. DETERMINE NONLINEAR TRANSFER. +! -------------------------------- + + IMPLICIT NONE + REAL, INTENT(IN) :: XI, XJ, XK, XIJ, XIK, XJK, XOI, XOJ, XOK + REAL :: RI, RJ, RK, OI, OJ, OK, SQIJK, SQIKJ, SQJKI + + RI=ABS(XI)+DEL1 + RJ=ABS(XJ)+DEL1 + RK=ABS(XK)+DEL1 + OI=XOI+DEL1 + OJ=XOJ+DEL1 + OK=XOK+DEL1 + SQIJK=SQRT(OI*OJ*RK/(OK*RI*RJ)) + SQIKJ=SQRT(OI*OK*RJ/(OJ*RI*RK)) + SQJKI=SQRT(OJ*OK*RI/(OI*RJ*RK)) + VPLUS_D=ZCONST*( (XIJ+RI*RJ)*SQIJK + (XIK+RI*RK)*SQIKJ & + & + (XJK+RJ*RK)*SQJKI ) + + END FUNCTION VPLUS_D +! ----------------------------------------------------------------- + + END SUBROUTINE SECONDHH + !/ ------------------------------------------------------------------- / + !/ + !> + !> @brief Determines skewness paramters in order to obtain + !> correction on altimeter wave height + !> + !> @details Evaluate deviations from gaussianity following the work + !> of Srokosz and Longuet-Higgins. For second order + !> corrections to surface elevation, the approach of + !> Zaharov has been used. + !> + !> @param[in] NKHF Extended number of frequencies. + !> @param[out] FAC0 2nd order coef correction. + !> @param[out] FAC1 2nd order coef correction. + !> @param[out] FAC2 2nd order coef correction. + !> @param[out] FAC3 2nd order coef correction. + !> + !> @author P. Janssen @date 29-Mar-2024 + !> + SUBROUTINE SKEWNESS(A) + +!-------------------------------------------------------------------- + +!*****SKEWNESS** COMPUTES PARAMETERS OF THE NEARLY-GAUSSIAN +! DISTRIBUTION OF OCEAN WAVES AT A FIXED GRID POINT. + +! P.JANSSEN JULY 1997 + +! PURPOSE +! ------- +! DETERMINES SKEWNESS PARAMETERS IN ORDER TO OBTAIN +! CORRECTION ON ALTIMETER WAVE HEIGHT. + +! INTERFACE +! --------- +! *CALL* *SKEWNESS(IU06,F1,NCOLL,XKAPPA1,DELH_ALT)* + + + +! METHOD +! ------ +! EVALUATE DEVIATIONS FROM GAUSSIANITY FOLLOWING THE WORK +! OF SROKOSZ AND LONGUET-HIGGINS. FOR SECOND ORDER +! CORRECTIONS TO SURFACE ELEVATION THE APPROACH OF +! ZAKHAROV HAS BEEN USED. + +! EXTERNALS +! --------- +! NONE + +! REFERENCES +! ---------- +! M.A. SROKOSZ, J.G.R.,91,995-1006(1986) +! V.E. ZAKHAROV, HAMILTONIAN APPROACH(1967) +!-------------------------------------------------------------------- + + + +!-------------------------------------------------------------------- +! *TH* REAL DIRECTIONS IN RADIANS. +USE CONSTANTS, ONLY: GRAV, TPI, TPIINV +USE W3GDATMD, ONLY: NK, NTH, XFR, SIG, DTH, ECOS, ESIN, NSEAL +USE W3PARALL, ONLY: INIT_GET_ISEA +USE W3ADATMD, ONLY: CG, SKEW, EMBIA1, EMBIA2 + + + IMPLICIT NONE + + REAL, INTENT(IN) :: A(NTH,NK,0:NSEAL) + + INTEGER :: NKHF + REAL(KIND=4), DIMENSION(:,:,:,:) , ALLOCATABLE:: FAC0,FAC1,FAC2,FAC3 + + INTEGER :: M, K, M1, K1, M2, K2, I, J + INTEGER :: MSTART, JSEA + + REAL(KIND=4) :: CONX, DELTA + REAL(KIND=4) :: FH, DELF, XK1 + REAL(KIND=4) :: XPI, XPJ, XPK, XN, XFAC, CO1 + REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: F2 + REAL(KIND=4), DIMENSION(0:3,0:2,0:2) :: XMU, XLAMBDA + REAL(KIND=4), DIMENSION(:) , ALLOCATABLE:: SIGHF, DFIMHF, FAK + +! ---------------------------------------------------------------------- + + NKHF=NK+13 ! same offset as in ECWAM + + ALLOCATE(FAC0(NTH,NTH,NKHF,NKHF)) + ALLOCATE(FAC1(NTH,NTH,NKHF,NKHF)) + ALLOCATE(FAC2(NTH,NTH,NKHF,NKHF)) + ALLOCATE(FAC3(NTH,NTH,NKHF,NKHF)) + + CALL SECONDHH(NKHF,FAC0,FAC1,FAC2,FAC3) + + ALLOCATE(F2(NTH,NKHF)) + ALLOCATE(SIGHF(NKHF), DFIMHF(NKHF), FAK(NKHF)) + +! 1. COMPUTATION OF FREQUENCY-DIRECTION INCREMENT +! ----------------------------------------------- + + MSTART = 1 + + +#ifdef W3_OMPG + !$OMP PARALLEL DO PRIVATE(JSEA) +#endif + DO JSEA=1, NSEAL + XMU(:,:,:) = 0.0 + DO K=1,NTH + DO M=1,NK + CONX = TPIINV / SIG(M) * CG(M,JSEA) + F2(K,M)=A(K,M,JSEA)/ CONX + END DO + END DO + + SIGHF(1) = SIG(1) + DO M=2,NKHF + SIGHF(M) = XFR*SIGHF(M-1) + ENDDO + + CO1 = 0.5*(XFR-1.)*DTH*TPIINV + DFIMHF(1) = CO1*SIGHF(1) ! this is DF*DTH + DO M=2,NKHF-1 + DFIMHF(M)=CO1*(SIGHF(M)+SIGHF(M-1)) + ENDDO + DFIMHF(NKHF)=CO1*SIGHF(NKHF-1) + + DO M=1,NKHF + FAK(M) = (SIGHF(M))**2/GRAV + ENDDO + +! Deals with the tail ... + DO M=NK+1,NKHF + FH=(SIGHF(NK)/SIGHF(M))**5 + DO K=1,NTH + F2(K,M)=F2(K,NK)*FH + ENDDO + ENDDO + +! 2. COMPUTATION OF THE SKEWNESS COEFFICIENTS +! -------------------------------------------- + + DO M1=MSTART,NKHF + DO M2=MSTART,NKHF + DO K1=1,NTH + DO K2=1,NTH + DELF = DFIMHF(M1)*DFIMHF(M2)*F2( K1,M1)*F2(K2,M2) + XMU(3,0,0) = XMU(3,0,0)+3.0*FAC0(K1,K2,M1,M2)*DELF + XMU(1,2,0) = XMU(1,2,0)+FAC1(K1,K2,M1,M2)*DELF + XMU(1,0,2) = XMU(1,0,2)+FAC2(K1,K2,M1,M2)*DELF + XMU(1,1,1) = XMU(1,1,1)+FAC3(K1,K2,M1,M2)*DELF + ENDDO + ENDDO + ENDDO + ENDDO + + DO K1=1,NTH + DO M1=MSTART,NKHF + XK1 = FAK(M1)**2 + DELF = DFIMHF(M1)*F2(K1,M1) + XMU(2,0,0) = XMU(2,0,0) + DELF + XMU(0,2,0) = XMU(0,2,0) + XK1*ECOS(K1)**2*DELF + XMU(0,0,2) = XMU(0,0,2) + XK1*ESIN(K1)**2*DELF + XMU(0,1,1) = XMU(0,1,1) + XK1*ECOS(K1)*ESIN(K1)*DELF + ENDDO + ENDDO + + +! 3. COMPUTATION OF THE NORMALISED SKEWNESS COEFFICIENTS +! ------------------------------------------------------ + + DO I=0,3 + XPI = 0.5*FLOAT(I) + DO J=0,2 + XPJ = 0.5*FLOAT(J) + DO K=0,2 + XPK = 0.5*FLOAT(K) + XN = XMU(2,0,0)**XPI*XMU(0,2,0)**XPJ*XMU(0,0,2)**XPK ! denom in Srokosz eq. 11 + IF (XN .NE. 0) THEN + XLAMBDA(I,J,K) = XMU(I,J,K)/XN + ELSE + XLAMBDA(I,J,K) = 0 + END IF + END DO + END DO + END DO + IF ( XMU(2,0,0) .GT. 1.E-7 ) THEN + SKEW(JSEA)=XLAMBDA(3,0,0) + DELTA = ( XLAMBDA(1,2,0) + XLAMBDA(1,0,2) & + - 2.0*XLAMBDA(0,1,1)*XLAMBDA(1,1,1) )/ & + (1.0 - XLAMBDA(0,1,1)**2) ! this is called gamma eq. 20 + EMBIA1(JSEA)=-0.125*DELTA ! EM Bias coefficient + EMBIA2(JSEA)=-0.125*XLAMBDA(3,0,0)/3.0 ! tracker bias (least squares only) + END IF + END DO ! end of loop on JSEA + ! +#ifdef W3_OMPG + !$OMP END PARALLEL DO +#endif + + DEALLOCATE(FAC0,FAC1,FAC2,FAC3) + DEALLOCATE(F2,SIGHF,DFIMHF,FAK) + + + END SUBROUTINE SKEWNESS + END MODULE W3IOGOMD diff --git a/model/src/w3odatmd.F90 b/model/src/w3odatmd.F90 index d268793fb..3a667ebbf 100644 --- a/model/src/w3odatmd.F90 +++ b/model/src/w3odatmd.F90 @@ -887,7 +887,7 @@ SUBROUTINE W3NOUT ( NDSERR, NDSTST ) ! ! 8) Spectrum parameters ! - NOGE(8) = 6 + NOGE(8) = 9 ! IDOUT( 8, 1) = 'Mean square slopes ' IDOUT( 8, 2) = 'Phillips tail const' @@ -895,6 +895,9 @@ SUBROUTINE W3NOUT ( NDSERR, NDSTST ) IDOUT( 8, 4) = 'Tail slope direction' IDOUT( 8, 5) = 'Goda peakedness parm' IDOUT( 8, 6) = 'kxky-peakdness ' + IDOUT( 8, 7) = 'Skewness ' + IDOUT( 8, 8) = 'EM bias(l120+l102)/8' + IDOUT( 8, 9) = 'Tracker bias:-l300/8' ! IDOUT( 8, 3) = 'Lx-Ly mean wvlength' ! IDOUT( 8, 4) = 'Surf grad correl XT' ! IDOUT( 8, 5) = 'Surf grad correl YT' diff --git a/model/src/w3ounfmetamd.F90 b/model/src/w3ounfmetamd.F90 index a4a58d079..87e606e56 100644 --- a/model/src/w3ounfmetamd.F90 +++ b/model/src/w3ounfmetamd.F90 @@ -3969,6 +3969,45 @@ SUBROUTINE DEFAULT_META() META(1)%VARNC='2D wavenumber peakedness' META(1)%VMIN = 0 META(1)%VMAX = 1600 + ! IFI=8, IFJ=7, SKW + META => GROUP(8)%FIELD(7)%META + META(1)%FSC = 0.00001 + META(1)%UNITS = '1' + META(1)%ENAME = '.skw' + META(1)%VARNM='skw' + META(1)%VARNL='skewness' + !META(1)%VARNS='sea_surface_wave_peakedness' + META(1)%VARNS='' + META(1)%VARNG='skewness of P(z,sx,sy=0)' + META(1)%VARNC='skewness of P(z,sx,sy=0)' + META(1)%VMIN = 0 + META(1)%VMAX = 1 + ! IFI=8, IFJ=8, EMB + META => GROUP(8)%FIELD(8)%META + META(1)%FSC = 0.00001 + META(1)%UNITS = '1' + META(1)%ENAME = '.emb' + META(1)%VARNM='emb' + META(1)%VARNL='EM-bias' + !META(1)%VARNS='sea_surface_wave_peakedness' + META(1)%VARNS='' + META(1)%VARNG='EM bias coefficient' + META(1)%VARNC='EM bias coefficient' + META(1)%VMIN = -1 + META(1)%VMAX = 1 + ! IFI=8, IFJ=7, SKW + META => GROUP(8)%FIELD(9)%META + META(1)%FSC = 0.00001 + META(1)%UNITS = '1' + META(1)%ENAME = '.emc' + META(1)%VARNM='emc' + META(1)%VARNL='trackerbias' + !META(1)%VARNS='sea_surface_wave_peakedness' + META(1)%VARNS='' + META(1)%VARNG='tracker bias coefficient' + META(1)%VARNC='tracker bias coefficient' + META(1)%VMIN = -1 + META(1)%VMAX = 1 ! ! !---------- GROUP 9 ---------------- ! diff --git a/model/src/ww3_ounf.F90 b/model/src/ww3_ounf.F90 index 02fd0d6f8..a2ff83e26 100644 --- a/model/src/ww3_ounf.F90 +++ b/model/src/ww3_ounf.F90 @@ -66,6 +66,7 @@ PROGRAM W3OUNF !/ 22-Mar-2021 : New coupling fields output ( version 7.12 ) !/ 02-Sep-2021 : Added coordinates attribute ( version 7.12 ) !/ 14-Feb-2023 : Added QKK output ( version 7.12 ) + !/ 03-Mar-2024 : Added SKEW & EMBIAS output ( version 7.xx ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -193,7 +194,7 @@ PROGRAM W3OUNF CFLTHMAX, CFLXYMAX, CFLKMAX, TAUICE, PHICE, & STMAXE, STMAXD, HMAXE, HCMAXE, HMAXD, HCMAXD,& P2SMS, EF, US3D, TH1M, STH1M, TH2M, STH2M, & - WN, USSP, WBT, WNMEAN, QKK + WN, USSP, WBT, WNMEAN, QKK, SKEW, EMBIA1, EMBIA2 USE W3ODATMD, ONLY: NDSO, NDSE, SCREEN, NOGRP, NGRPP, IDOUT, & UNDEF, FLOGRD, FNMPRE, NOSWLL, NOGE ! @@ -1960,6 +1961,18 @@ SUBROUTINE W3EXNC ( NX, NY, IX1, IXN, IY1, IYN, NSEA, & ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 6 ) THEN CALL S2GRID(QKK, X1) ! + ! surface elevation skewness lambda_3,0,0 + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 7 ) THEN + CALL S2GRID(SKEW, X1) + ! + ! em bias param 1 + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 8 ) THEN + CALL S2GRID(EMBIA1, X1) + ! + ! em bias param 2 + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 9 ) THEN + CALL S2GRID(EMBIA2, X1) + ! ! Dynamic time step ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 1 ) THEN DO ISEA=1, NSEA diff --git a/model/src/ww3_outf.F90 b/model/src/ww3_outf.F90 index e955f00ad..590518037 100644 --- a/model/src/ww3_outf.F90 +++ b/model/src/ww3_outf.F90 @@ -160,6 +160,7 @@ PROGRAM W3OUTF PHS, PTP, PLP, PDIR, PSI, PWS, PWST, PNR, & PTM1, PT1, PT2, PEP, TAUOCX, TAUOCY, & PTHP0, PQP, PSW, PPE, PGW, QP, QKK, & + SKEW, EMBIA1, EMBIA2, & TAUOX, TAUOY, TAUWIX,BHD, & TAUWIY, PHIAW, PHIOC, TUSX, TUSY, PRMS, TPMS,& USSX, USSY, MSSX, MSSY, MSCX, MSCY, CHARN, & @@ -2216,6 +2217,39 @@ SUBROUTINE W3EXGO ( NX, NY, NSEA ) CALL W3S2XY ( NSEA, NSEA, NX+1, NY, QKK, MAPSF, X1 ) ENDIF ! + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 7 ) THEN + FLONE = .TRUE. + FSC = 0.01 + UNITS = '1' + ENAME = '.skw' + IF ( ITYPE .EQ. 4 ) THEN + XS1 = SKEW + ELSE + CALL W3S2XY ( NSEA, NSEA, NX+1, NY, SKEW, MAPSF, X1 ) + ENDIF + ! + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 8 ) THEN + FLONE = .TRUE. + FSC = 0.0001 + UNITS = '1' + ENAME = '.emb' + IF ( ITYPE .EQ. 4 ) THEN + XS1 = EMBIA1 + ELSE + CALL W3S2XY ( NSEA, NSEA, NX+1, NY, EMBIA1, MAPSF, X1 ) + ENDIF + ! + ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 9 ) THEN + FLONE = .TRUE. + FSC = 0.0001 + UNITS = '1' + ENAME = '.emc' + IF ( ITYPE .EQ. 4 ) THEN + XS1 = EMBIA2 + ELSE + CALL W3S2XY ( NSEA, NSEA, NX+1, NY, EMBIA2, MAPSF, X1 ) + ENDIF + ! ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 1 ) THEN FLONE = .TRUE. FSC = 0.1 diff --git a/model/tools/bash/ww3_ounf_inp2nml.sh b/model/tools/bash/ww3_ounf_inp2nml.sh index e9a34d7b8..1edd05562 100755 --- a/model/tools/bash/ww3_ounf_inp2nml.sh +++ b/model/tools/bash/ww3_ounf_inp2nml.sh @@ -184,7 +184,7 @@ cat >> $nmlfile << EOF ! UST CHA CGE FAW TAW TWA WCC WCF WCH WCM FWS ! SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC USP TOC ! ABR UBR BED FBB TBB -! MSS MSC WL02 AXT AYT AXY +! MSS MSC MSD MCD QP QKK SKW EMB EMC ! DTD FC CFX CFD CFK ! U1 U2 ! diff --git a/model/tools/bash/ww3_shel_inp2nml.sh b/model/tools/bash/ww3_shel_inp2nml.sh index 619002aa8..8ea336e13 100755 --- a/model/tools/bash/ww3_shel_inp2nml.sh +++ b/model/tools/bash/ww3_shel_inp2nml.sh @@ -970,6 +970,9 @@ cat >> $nmlfile << EOF ! F F 8 4 MSCD MCD Tail slope direction ! F F 8 5 QP QP Goda peakedness parameter ! F F 8 6 QKK QKK Wavenumber peakedness +! F F 8 7 SKEW SKW Skewness of elevation for zero slopes +! F F 8 8 EMBIA1 EMB Mean sea level at zero slopes / Hs +! F F 8 9 EMBIA2 EMC Tracker bias for LRM least square altimetry ! ------------------------------------------------- ! 9 Numerical diagnostics ! ------------------------------------------------- diff --git a/regtests/ww3_ts1/input/ww3_ounf.inp b/regtests/ww3_ts1/input/ww3_ounf.inp index 52a2dd2c6..0b89ef01e 100644 --- a/regtests/ww3_ts1/input/ww3_ounf.inp +++ b/regtests/ww3_ts1/input/ww3_ounf.inp @@ -11,7 +11,7 @@ $ file for a full documentation of field output options. Namelist type $ selection is used here (for alternative F/T flags, see ww3_shel.inp). $ N - DPT WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT + DPT WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT SKW EMB EMC $ $--------------------------------------------------------------------- $ $ NetCDF version [3,4] and variable type 4 [2 = SHORT, 3 = it depends , 4 = REAL] diff --git a/regtests/ww3_ts1/input/ww3_ounf.nml b/regtests/ww3_ts1/input/ww3_ounf.nml index fb0f02d3e..f9992f0ce 100644 --- a/regtests/ww3_ts1/input/ww3_ounf.nml +++ b/regtests/ww3_ts1/input/ww3_ounf.nml @@ -9,7 +9,7 @@ FIELD%TIMESTART = '19680101 120000' FIELD%TIMESTRIDE = '10' FIELD%TIMECOUNT = '8000' - FIELD%LIST = 'DPT WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT' + FIELD%LIST = 'DPT WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT SKW EMB EMC' FIELD%PARTITION = '0 1 2' FIELD%TYPE = 4 / diff --git a/regtests/ww3_ts1/input/ww3_shel.inp b/regtests/ww3_ts1/input/ww3_shel.inp index fca96fb7e..5171bd9ab 100644 --- a/regtests/ww3_ts1/input/ww3_shel.inp +++ b/regtests/ww3_ts1/input/ww3_shel.inp @@ -19,7 +19,7 @@ $ $ N $ -DPT WND MSS MSD ICE HS MSS FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT +DPT WND MSS MSD ICE HS MSS FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT SKW EMB EMC $ 19680606 000000 60 19680618 000000 0.0 0.0 'The_point' diff --git a/regtests/ww3_ts1/input/ww3_shel.nml b/regtests/ww3_ts1/input/ww3_shel.nml index b4837d6e1..1ecf48c51 100644 --- a/regtests/ww3_ts1/input/ww3_shel.nml +++ b/regtests/ww3_ts1/input/ww3_shel.nml @@ -21,7 +21,7 @@ ! Define the output types point parameters via OUTPUT_TYPE_NML namelist ! -------------------------------------------------------------------- ! &OUTPUT_TYPE_NML - TYPE%FIELD%LIST = 'DPT WND MSS MSD ICE HS MSS FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT' + TYPE%FIELD%LIST = 'DPT WND MSS MSD ICE HS MSS FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT SKW EMB EMC' TYPE%POINT%FILE = '../input/points.list' /