!#include "w3macros.h"
#include PROJECT_HEADER
#include GLOBAL_DEFS
!/ ------------------------------------------------------------------- /
MODULE W3UPDTMD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 13-Jan-2016 |
!/ +-----------------------------------+
!/
!/ 21-Jan-2000 : Origination. ( version 2.00 )
!/ 24-Jan-2001 : Flat grid version. ( version 2.06 )
!/ 02-Apr-2001 : Adding sub-grid obstacles. ( version 2.10 )
!/ 18-May-2001 : Clean up and bug fixes. ( version 2.11 )
!/ 11-Jan-2002 : Sub-grid ice. ( version 2.15 )
!/ 30-Apr-2002 : Water level fixes. ( version 2.20 )
!/ 13-Nov-2002 : Add stress vector. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid wind correction. ( version 3.02 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 15-Jul-2005 : Adding MAPST2. ( version 3.07 )
!/ 07-Sep-2005 : Upgrading W3UBPT. ( version 3.08 )
!/ 04-Jul-2006 : Consolidate stress arrays. ( version 3.09 )
!/ 11-Jan-2007 : Clean-up W3UTRN boundary points. ( version 3.10 )
!/ 11-May-2007 : Adding NTPROC/NAPROC separation. ( version 3.11 )
!/ 29-May-2009 : Preparing distribution version. ( version 3.14 )
!/ 30-Oct-2009 : Implement run-time grid selection. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 17-Aug-2010 : ABPI0-N(:,0) init. bug fix. ( version 3.14 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/ 05-Apr-2011 : Place holder for XGR in UNGTYPE ( version 4.04 )
!/ (A. Roland/F. Ardhuin)
!/ 13-Mar-2012 : Add initialization of UST on re- ( version 4.07 )
!/ activation of grid point.
!/ 06-Jun-2012 : Porting bugfixes from 3.14 to 4.07 ( version 4.07 )
!/ 12-Jun-2012 : Add /RTD option or rotated grid option.
!/ (Jian-Guo Li) ( version 4.07 )
!/ 26-Sep-2012 : Adding update from tidal analysis ( version 4.08 )
!/ (F. Ardhuin)
!/ 16-Sep-2013 : Add Arctic part for SMC grid. ( version 4.11 )
!/ 11-Nov-2013 : SMC and rotated grid incorporated in the main
!/ trunk ( version 4.13 )
!/ 13-Nov-2013 : Moved reflection from ww3_grid.ftn ( version 4.13 )
!/ 27-May-2014 : Ading OMPG parallelizations dir, ( version 5.02 )
!/ 08-May-2014 : Implement tripolar grid for first order propagation
!/ scheme ( version 5.03 )
!/ (W. E. Rogers, NRL)
!/ 27-Aug-2015 : New function to update variables ( version 5.08 )
!/ ICEF and ICEDMAX at the first time step
!/ and add ICEH initialization in W3UICE.
!/ 13-Jan-2016 : Changed initial value of ICEDMAX ( version 5.08 )
!/
!/ Copyright 2009-2014 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 :
!
! Bundles all input updating routines for WAVEWATCH III.
!
! 2. Variables and types :
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3UCUR Subr. Public Update current fields.
! W3UWND Subr. Public Update wind fields.
! W3UINI Subr. Public Update initial conditions.
! W3UBPT Subr. Public Update boundary conditions.
! W3UICE Subr. Public Update ice concentrations.
! W3ULEV Subr. Public Update water levels.
! W3UTRN Subr. Public Update cell boundary transparancies.
! W3DZXY Subr. Public Calculate derivatives of a field.
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! DSEC21 Func. W3TIMEMD Difference in time.
! STRACE Subr. W3SERVMD Subroutine tracing.
! EXTCDE Subr. W3SERVMD Exit program with error code.
! PRTBLK Subr. W3ARRYMD Print plot output.
! PRT2DS Subr. W3ARRYMD Print plot output.
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! 6. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/DIST Id.
!
! !/OMPG OpenMP compiler directives.
!
! !/CRT0 No current interpolation.
! !/CRT1 Linear current interpolation.
! !/CRT2 Quasi-quadratic current interpolation.
!
! !/WNT0 No wind interpolation.
! !/WNT1 Linear wind interpolation.
! !/WNT2 Energy conservation in wind interpolation.
!
! !/RWND Use wind speeds relative to currents.
!
! !/STAB2 Calculate effective wind speed factor for stability
! to be used with !/ST2.
!
! !/S Enable subroutine tracing.
! !/Tn Test output
!
! 7. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
USE W3ODATMD, ONLY: NDSE, NDST, NAPROC, IAPROC, NAPERR
!/S USE W3SERVMD, ONLY : STRACE
USE W3TIMEMD, ONLY: DSEC21
!/
!/ ------------------------------------------------------------------- /
!/
CONTAINS
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UCUR ( FLFRST )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 15-Dec-2004 |
!/ +-----------------------------------+
!/
!/ 09-Dec-1996 : Final FORTRAN 77 ( version 1.18 )
!/ 20-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 27-Aug-2015 : Rename DT0,DTT by DT0T,DT0N ( version 5.10 )
!/ 23-Mar-2016 : SMC grid Arctic part adjustment. ( version 5.18 )
!/
! 1. Purpose :
!
! Interpolate the current field to the present time.
!
! 2. Method :
!
! Linear interpolation of speed and direction, with optionally
! a correction to get approximate quadratic interpolation of speed
! only.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! FLFRST Log. I Flag for first pass through routine.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Only currents at sea points are considered.
! - Time ranges checked in W3WAVE.
! - Currents are stored by components to save on the use of
! SIN and COS functions. The actual interpolations, however
! are by absolute value and direction.
!
! 8. Structure :
!
! --------------------------------------
! 1. Prepare auxiliary arrays.
! 2. Calculate interpolation factors.
! 3. Get actual winds.
! --------------------------------------
!
! 9. Switches :
!
! !/CRT0 No current interpolation.
! !/CRT1 Linear current interpolation.
! !/CRT2 Quasi-quadratic current interpolation.
!
! !/S Enable subroutine tracing.
! !/T Test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
!/ARC USE W3GDATMD, ONLY: NARC, NGLO, ANGARC
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: CX, CY, CA0, CAI, CD0, CDI
USE W3IDATMD, ONLY: TC0, CX0, CY0, TCN, CXN, CYN
!/TIDE USE W3GDATMD, ONLY: YGRD
!/TIDE USE W3TIMEMD
!/TIDE USE W3IDATMD, ONLY: FLCURTIDE, CXTIDE, CYTIDE, NTIDE
!/TIDE USE W3TIDEMD
!
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
LOGICAL, INTENT(IN) :: FLFRST
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: ISEA, IX, IY
!/S INTEGER, SAVE :: IENT = 0
REAL :: D0, DN, DD, DT0N, DT0T, RD, CABS, CDIR
!/CRT2 REAL :: RD2, CI2
!/TIDE INTEGER :: J,K
!/TIDE INTEGER(KIND=4) :: TIDE_KD0, INT24, INTDYS ! "Gregorian day constant"
!/TIDE REAL :: WCURTIDEX, WCURTIDEY, TIDE_ARGX, TIDE_ARGY
!/TIDE REAL(KIND=8) :: d1,h,TIDE_HOUR,HH,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau
!/TIDE REAL :: FX(44),UX(44),VX(44)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UCUR')
!
! 1. Prepare auxiliary arrays
!
IF ( FLFRST ) THEN
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
CA0(ISEA) = SQRT ( CX0(IX,IY)**2 + CY0(IX,IY)**2 )
CAI(ISEA) = SQRT ( CXN(IX,IY)**2 + CYN(IX,IY)**2 )
IF ( CA0(ISEA) .GT. 1.E-7) THEN
D0 = MOD ( TPI+ATAN2(CY0(IX,IY),CX0(IX,IY)) , TPI )
ELSE
D0 = 0
END IF
IF ( CAI(ISEA) .GT. 1.E-7) THEN
DN = MOD ( TPI+ATAN2(CYN(IX,IY),CXN(IX,IY)) , TPI )
ELSE
DN = D0
END IF
IF ( CA0(ISEA) .GT. 1.E-7) THEN
CD0(ISEA) = D0
ELSE
CD0(ISEA) = DN
END IF
DD = DN - CD0(ISEA)
IF (ABS(DD).GT.PI) DD = DD - TPI*SIGN(1.,DD)
CDI(ISEA) = DD
CAI(ISEA) = CAI(ISEA) - CA0(ISEA)
END DO
END IF
!
! 2. Calculate interpolation factor
!
DT0N = DSEC21 ( TC0, TCN )
DT0T = DSEC21 ( TC0, TIME )
!
!/CRT0 RD = 0.
!/CRT1 RD = DT0T / MAX ( 1.E-7 , DT0N )
!/CRT2 RD = DT0T / MAX ( 1.E-7 , DT0N )
!/CRT2 RD2 = 1. - RD
!/OASIS RD = 1.
!
!/T WRITE (NDST,9000) DT0N, DT0T, RD
!/TIDE IF (FLCURTIDE) THEN
!/TIDE! WRITE(6,*) 'TIME CUR:',TIME, '##',TC0, '##',TCN
!/TIDE TIDE_HOUR = TIME2HOURS(TIME)
!/TIDE!
!/TIDE!* THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
!/TIDE!* AT THE MID POINT OF THE ANALYSIS PERIOD.
!/TIDE d1=TIDE_HOUR/24.d0
!/TIDE TIDE_KD0= 2415020
!/TIDE d1=d1-dfloat(TIDE_kd0)-0.5d0
!/TIDE call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
!/TIDE INT24=24
!/TIDE INTDYS=int((TIDE_HOUR+0.00001)/INT24)
!/TIDE HH=TIDE_HOUR-dfloat(INTDYS*INT24)
!/TIDE TAU=HH/24.D0+H-S
!/TIDE END IF
!/TIDE!
!/TIDE! ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
!/TIDE! TING THE LUNAR TIME TAU.
!/TIDE!
!
! 3. Actual currents for all grid points
!
DO ISEA=1, NSEA
!/TIDE IF (FLCURTIDE) THEN ! could move IF test outside of ISEA loop ...
!/TIDE! VUF should only be updated in latitude changes significantly ...
!/TIDE IX = MAPSF(ISEA,1)
!/TIDE IY = MAPSF(ISEA,2)
!/TIDE CALL SETVUF_FAST(h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau,YGRD(IY,IX),FX,UX,VX)
!/TIDE WCURTIDEX = CXTIDE(IX,IY,1,1)
!/TIDE WCURTIDEY = CYTIDE(IX,IY,1,1)
!/TIDE
!/TIDE DO J=2,TIDE_MF
!/TIDE TIDE_ARGX=(VX(J)+UX(J))*twpi-CXTIDE(IX,IY,J,2)*DERA
!/TIDE TIDE_ARGY=(VX(J)+UX(J))*twpi-CYTIDE(IX,IY,J,2)*DERA
!/TIDE WCURTIDEX = WCURTIDEX+FX(J)*CXTIDE(IX,IY,J,1)*COS(TIDE_ARGX)
!/TIDE WCURTIDEY = WCURTIDEY+FX(J)*CYTIDE(IX,IY,J,1)*COS(TIDE_ARGY)
!/TIDE END DO
!/TIDE
!/TIDET !Verification
!/TIDET IF (ISEA.EQ.1) THEN
!/TIDET
!/TIDET TIDE_AMPC(1:NTIDE,1)=CXTIDE(IX,IY,1:NTIDE,1)
!/TIDET TIDE_PHG(1:NTIDE,1 )=CXTIDE(IX,IY,1:NTIDE,2)
!/TIDET TIDE_AMPC(1:NTIDE,2)=CYTIDE(IX,IY,1:NTIDE,1)
!/TIDET TIDE_PHG(1:NTIDE,2) =CYTIDE(IX,IY,1:NTIDE,2)
!/TIDET
!/TIDET WRITE(993,'(A,F20.2,13F8.3)') 'TEST ISEA 0:', &
!/TIDET d1,H,S,TAU,pp,s,p,enp,dh,dpp,ds,dp,dnp,YGRD(IY,IX)
!/TIDET
!/TIDET DO J=1,TIDE_MF
!/TIDET WRITE(993,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',IX,J,TIME,TIDE_HOUR, &
!/TIDET FX(J),UX(J),VX(J),TIDE_INDEX2(J),TIDECON_ALLNAMES(TIDE_INDEX2(J))
!/TIDET END DO
!/TIDET DO K=1,2
!/TIDET DO J=1,TIDE_MF
!/TIDET WRITE(993,'(A,5I9,F12.0,5F8.3)') 'TEST ISEA 2:',IX,K,J,TIME,TIDE_HOUR, &
!/TIDET FX(J),UX(J),VX(J),TIDE_AMPC(J,K),TIDE_PHG(J,K)
!/TIDET END DO
!/TIDET END DO
!/TIDET
!/TIDET WRITE(993,'(A,2F8.4,A,2F8.4)') '#:',CX0(IX,IY),CY0(IX,IY),'##',WCURTIDEX,WCURTIDEY
!/TIDET CLOSE(993)
!/TIDET END IF
!/TIDET ! End of verification
!/TIDE CX(ISEA) = WCURTIDEX
!/TIDE CY(ISEA) = WCURTIDEY
!/TIDE ELSE
CABS = CA0(ISEA) + RD * CAI(ISEA)
!/CRT2 CI2 = SQRT ( RD2 * CA0(ISEA)**2 + &
!/CRT2 RD *(CA0(ISEA)+CAI(ISEA))**2 )
!/CRT2 CABS = CABS * MIN( 1.25 , CI2/MAX(1.E-7,CABS) )
CDIR = CD0(ISEA) + RD * CDI(ISEA)
!/ARC !Li Rotate curreent direction by ANGARC for Arctic part cells. JGLi23Mar2016
!/ARC IF( ISEA .GT. NGLO ) THEN
!/ARC DN = CDIR + ANGARC( ISEA - NGLO )*DERA
!/ARC CDIR = MOD ( TPI + DN, TPI )
!/ARC ENDIF
CX(ISEA) = CABS * COS(CDIR)
CY(ISEA) = CABS * SIN(CDIR)
!/TIDE ! IF (ISEA.EQ.1) WRITE(6,'(A,4F8.4,A,4F8.4)') 'CUR#:',RD,CA0(ISEA),CAI(ISEA),CABS,'##', &
!/TIDE ! CX(ISEA), CY(ISEA),WCURTIDEX, WCURTIDEY
!/TIDE END IF
!
END DO
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3UCUR : DT0N, DT0T, RD :',2F8.1,F6.3)
!/
!/ End of W3UCUR ----------------------------------------------------- /
!/
END SUBROUTINE W3UCUR
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UWND ( FLFRST, VGX, VGY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 27-May-2014 |
!/ +-----------------------------------+
!/
!/ 03-Dec-1998 : Final FORTRAN 77 ( version 1.18 )
!/ 20-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 13-Nov-2002 : Add stress vector. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid wind correction. ( version 3.02 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 04-Jul-2006 : Consolidate stress arrays. ( version 3.09 )
!/ 16-Sep-2013 : Rotating wind for Arctic part. ( version 4.11 )
!/ 27-May-2014 : Adding OMPG parallelizations dir. ( version 5.02 )
!/ 27-Aug-2015 : Rename DT0,DTT by DT0T,DT0N ( version 5.10 )
!/
! 1. Purpose :
!
! Interpolate wind fields to the given time.
!
! 2. Method :
!
! Linear interpolation of wind speed and direction, with a simple
! correction to obtain quasi-conservation of energy.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! FLFRST Log. I Flag for first pass through routine.
! VGX/Y Real I Grid velocity (!/MGW)
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Only winds over sea points are considered.
! - Time ranges checked in W3WAVE.
!
! 8. Structure :
!
! --------------------------------------
! 1. Prepare auxiliary arrays.
! 2. Calculate interpolation factors
! 3. Get actual winds
! 4. Correct for currents
! 5. Convert to stresses
! 6. Stability correction
! --------------------------------------
!
! 9. Switches :
!
! !/OMPG OpenMP compiler directives.
!
! !/WNT0 No wind interpolation.
! !/WNT1 Linear wind interpolation.
! !/WNT2 Energy conservation in wind interpolation.
!
! !/RWND Use wind speeds relative to currents.
! !/MGW Moving grid wind correction.
!
! !/STAB2 Calculate effective wind speed factor for stability
! to be used with !/ST2.
!
! !/S Enable subroutine tracing.
! !/T Test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
!/WCOR USE W3GDATMD, ONLY: WWCOR
!/RWND USE W3GDATMD, ONLY: RWINDC
!/ST2 USE W3GDATMD, ONLY: ZWIND, OFSTAB, FFNG, FFPS, CCNG, CCPS, SHSTAB
!/ARC USE W3GDATMD, ONLY: NARC, NGLO, ANGARC
USE W3WDATMD, ONLY: TIME, ASF
USE W3ADATMD, ONLY: DW, CX, CY, UA, UD, U10, U10D, AS, &
UA0, UAI, UD0, UDI, AS0, ASI
USE W3IDATMD, ONLY: TW0, WX0, WY0, DT0, TWN, WXN, WYN, DTN, FLCUR
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: VGX, VGY
LOGICAL, INTENT(IN) :: FLFRST
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: ISEA, IX, IY
!/S INTEGER, SAVE :: IENT = 0
REAL :: D0, DN, DD, DT0N, DT0T, RD, UI2, &
UXR, UYR
!/WNT2 REAL :: RD2
!/STAB2 REAL :: STAB0, STAB, THARG1, THARG2, COR1, COR2
REAL :: UDARC
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UWND')
!
! 1. Prepare auxiliary arrays
!
IF ( FLFRST ) THEN
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
UA0(ISEA) = SQRT ( WX0(IX,IY)**2 + WY0(IX,IY)**2 )
UAI(ISEA) = SQRT ( WXN(IX,IY)**2 + WYN(IX,IY)**2 )
IF ( UA0(ISEA) .GT. 1.E-7) THEN
D0 = MOD ( TPI+ATAN2(WY0(IX,IY),WX0(IX,IY)) , TPI )
ELSE
D0 = 0
END IF
IF ( UAI(ISEA) .GT. 1.E-7) THEN
DN = MOD ( TPI+ATAN2(WYN(IX,IY),WXN(IX,IY)) , TPI )
ELSE
DN = D0
END IF
IF ( UA0(ISEA) .GT. 1.E-7) THEN
UD0(ISEA) = D0
ELSE
UD0(ISEA) = DN
END IF
DD = DN - UD0(ISEA)
IF (ABS(DD).GT.PI) DD = DD - TPI*SIGN(1.,DD)
UDI(ISEA) = DD
UAI(ISEA) = UAI(ISEA) - UA0(ISEA)
AS0(ISEA) = DT0(IX,IY)
ASI(ISEA) = DTN(IX,IY) - DT0(IX,IY)
END DO
END IF
!
! 2. Calculate interpolation factor
!
DT0N = DSEC21 ( TW0, TWN )
DT0T = DSEC21 ( TW0, TIME )
!
!/WNT0 RD = 0.
!/WNT1 RD = DT0T / MAX ( 1.E-7 , DT0N )
!/WNT2 RD = DT0T / MAX ( 1.E-7 , DT0N )
!/WNT2 RD2 = 1. - RD
!/OASIS RD = 1.
!
!/T WRITE (NDST,9000) DT0N, DT0T, RD
!
! 3. Actual wind for all grid points
!
!/OMPG/!$OMP PARALLEL DO PRIVATE (ISEA,UI2,UXR,UYR,UDARC)
!
DO ISEA=1, NSEA
!
! jcw call to stop wind interpolation.
#if !defined AIR_WAVES
UA(ISEA) = UA0(ISEA) + RD * UAI(ISEA)
!/WNT2 UI2 = SQRT ( RD2 * UA0(ISEA)**2 + &
!/WNT2 RD *(UA0(ISEA)+UAI(ISEA))**2 )
!/WNT2 UA(ISEA) = UA(ISEA) * MIN(1.25,UI2/MAX(1.E-7,UA(ISEA)))
UD(ISEA) = UD0(ISEA) + RD * UDI(ISEA)
#endif
!/MGW UXR = UA(ISEA)*COS(UD(ISEA)) + VGX
!/MGW UYR = UA(ISEA)*SIN(UD(ISEA)) + VGY
!/MGW UA(ISEA) = MAX ( 0.001 , SQRT(UXR**2+UYR**2) )
!/MGW UD(ISEA) = MOD ( TPI+ATAN2(UYR,UXR) , TPI )
!/ARC !Li Rotate wind direction by ANGARC for Arctic part cells.
!/ARC IF( ISEA .GT. NGLO ) THEN
!/ARC UDARC = UD(ISEA) + ANGARC( ISEA - NGLO )*DERA
!/ARC UD(ISEA) = MOD ( TPI + UDARC, TPI )
!/ARC ENDIF
!
AS(ISEA) = AS0(ISEA) + RD * ASI(ISEA)
! IF (UA(ISEA).NE.UA(ISEA)) WRITE(6,*) 'BUG WIND:',ISEA,UA(ISEA),MAPSF(ISEA,1), MAPSF(ISEA,2),UA0(ISEA),RD,UAI(ISEA)
! IF (UD(ISEA).NE.UD(ISEA)) WRITE(6,*) 'BUG WIN2:',ISEA,UD(ISEA),MAPSF(ISEA,1), MAPSF(ISEA,2)
!
END DO
!
!/OMPG/!$OMP END PARALLEL DO
!
! 3.b Bias correction ( !/WCOR )
!/WCOR WHERE ( UA .GE. WWCOR(1) ) UA = UA+(UA-WWCOR(1))*WWCOR(2)
!
! 4. Correct for currents and grid motion
!
!/RWND IF ( FLCUR ) THEN
!
!/RWND DO ISEA=1, NSEA
!/RWND UXR = UA(ISEA)*COS(UD(ISEA)) - RWINDC*CX(ISEA)
!/RWND UYR = UA(ISEA)*SIN(UD(ISEA)) - RWINDC*CY(ISEA)
!/RWND U10 (ISEA) = MAX ( 0.001 , SQRT(UXR**2+UYR**2) )
!/RWND U10D(ISEA) = MOD ( TPI+ATAN2(UYR,UXR) , TPI )
!/RWND END DO
!
!/RWND ELSE
!
!/OMPG/!$OMP PARALLEL DO PRIVATE (ISEA)
!
DO ISEA=1, NSEA
U10 (ISEA) = MAX ( UA(ISEA) , 0.001 )
U10D(ISEA) = UD(ISEA)
END DO
!
!/OMPG/!$OMP END PARALLEL DO
!
!/RWND END IF
!
! 5. Stability correction ( !/STAB2 )
! Original settings :
!
! SHSTAB = 1.4
! OFSTAB = -0.01
! CCNG = -0.1
! CCPS = 0.1
! FFNG = -150.
! FFPS = 150.
!
!/STAB2 STAB0 = ZWIND * GRAV / 273.
!
!/STAB2 DO ISEA=1, NSEA
!/STAB2 STAB = STAB0 * AS(ISEA) / MAX(5.,U10(ISEA))**2
!/STAB2 STAB = MAX ( -1. , MIN ( 1. , STAB ) )
!
!/STAB2 THARG1 = MAX ( 0. , FFNG*(STAB-OFSTAB))
!/STAB2 THARG2 = MAX ( 0. , FFPS*(STAB-OFSTAB))
!/STAB2 COR1 = CCNG * TANH(THARG1)
!/STAB2 COR2 = CCPS * TANH(THARG2)
!
!/STAB2 ASF(ISEA) = SQRT ( (1.+COR1+COR2)/SHSTAB )
!/STAB2 U10(ISEA) = U10(ISEA) / ASF(ISEA)
!/STAB2 END DO
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3UWND : DT0N, DT0T, RD :',2F8.1,F6.3)
!/
!/ End of W3UWND ----------------------------------------------------- /
!/
END SUBROUTINE W3UWND
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UINI ( A )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 05-Apr-2011 |
!/ +-----------------------------------+
!/
!/ 19-Oct-1998 : Final FORTRAN 77 ( version 1.18 )
!/ 20-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 24-Jan-2001 : Flat grid version. ( version 2.06 )
!/ 18-May-2001 : Fix CG declaration. ( version 2.11 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 30-Oct-2009 : Implement run-time grid selection. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/
! 1. Purpose :
!
! Initialize the wave field with fetch-limited spectra before the
! actual calculation start. (Named as an update routine due to
! placement in code.)
!
! 2. Method :
!
! Fetch-limited JONSWAP spectra with a cosine^2 directional
! distribution and a mean direction taken from the wind.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! A R.A. O Action density spectra.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Wind speeds filtered by U10MIN and U10MAX (DATA statements)
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/DIST Id.
!
! !/S Enable subroutine tracing.
! !/T General test output.
! !/T1 Parameters at grid points.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF, &
NK, NTH, TH, SIG, DTH, DSIP, UNGTYPE, &
RLGTYPE, CLGTYPE, GTYPE, FLAGLL, &
HPFAC, HQFAC
USE W3ADATMD, ONLY: U10, U10D, CG
!/T USE W3ARRYMD, ONLY : PRTBLK
!
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(OUT) :: A(NTH,NK,0:NSEAL)
!/
!/ ------------------------------------------------------------------- /
!/ Local variables
!/
INTEGER :: IX, IY, ISEA, JSEA, IK, ITH
!/S INTEGER, SAVE :: IENT = 0
!/T INTEGER :: IX0, IXN, MAPOUT(NX,NY)
!/T INTEGER :: NXP = 60
REAL :: ALFA(NSEAL), FP(NSEAL), YLN(NSEAL), &
AA, BB, CC
REAL :: XGR, U10C, U10DIR, XSTAR, FSTAR, &
GAMMA, FR, D1(NTH), D1INT, F1, F2
REAL :: ETOT, E1I
REAL :: U10MIN = 1.
REAL :: U10MAX = 20.
!/T REAL :: HSIG(NX,NY)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UINI')
!
!
! Pre-process JONSWAP data for all grid points ----------------------- *
!
!/T1 WRITE (NDST,9010)
!
! this is not clear what is going on betwen w3init and this ...
A(:,:,:)=0
DO JSEA=1, NSEAL
!/DIST ISEA = IAPROC + (JSEA-1)*NAPROC
!/SHRD ISEA = JSEA
!
IF (GTYPE.EQ.UNGTYPE) THEN
XGR=1. ! to be fixed later
ELSE
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
XGR = 0.5 * SQRT(HPFAC(IY,IX)**2+HQFAC(IY,IX)**2)
END IF
IF ( FLAGLL ) THEN
XGR = XGR * RADIUS * DERA
END IF
!
U10C = MAX ( MIN(U10(ISEA),U10MAX) , U10MIN )
!
XSTAR = GRAV * XGR / U10C**2
FSTAR = 3.5 / XSTAR**(0.33)
GAMMA = MAX ( 1. , 7.0 / XSTAR**(0.143) )
!
ALFA(JSEA) = 0.076 / XSTAR**(0.22)
FP (JSEA) = FSTAR * GRAV / U10C
YLN (JSEA) = LOG ( GAMMA )
!
!/T1 WRITE (NDST,9011) ISEA, U10C, XSTAR, &
!/T1 ALFA(JSEA), FP(JSEA), GAMMA
!
END DO
!
! 1-D spectrum at location ITH = NTH --------------------------------- *
!
DO IK=1, NK
FR = SIG(IK) * TPIINV
DO JSEA=1, NSEAL
!
!/ ----- INLINED EJ5P (REDUCED) -------------------------------------- /
!
AA = ALFA(JSEA) * 0.06175/FR**5
BB = MAX( -50. , -1.25*(FP(JSEA)/FR)**4 )
CC = MAX( -50. , -0.5*((FR-FP(JSEA))/(0.07*FP(JSEA)))**2 )
A(NTH,IK,JSEA) &
= AA * EXP(BB + EXP(CC) * YLN(JSEA))
!
!/ ----- INLINED EJ5P (END) ------------------------------------------ /
!
END DO
END DO
!
! Apply directional distribution ------------------------------------- *
!
DO JSEA=1, NSEAL
!/DIST ISEA = IAPROC + (JSEA-1)*NAPROC
!/SHRD ISEA = JSEA
!
U10DIR = U10D(ISEA)
D1INT = 0.
!
DO ITH=1, NTH
D1(ITH) = ( MAX ( 0. , COS(TH(ITH)-U10DIR) ) )**2
D1INT = D1INT + D1(ITH)
END DO
!
D1INT = D1INT * DTH
F1 = TPIINV / D1INT
!
DO IK=1, NK
F2 = F1 * A(NTH,IK,JSEA) * CG(IK,ISEA) / SIG(IK)
DO ITH=1, NTH
A(ITH,IK,JSEA) = F2 * D1(ITH)
END DO
END DO
!
END DO
!
! Test output -------------------------------------------------------- *
!
!/T HSIG = 0.
!/T MAPOUT = 0
!
!/T DO ISEA=IAPROC, NSEA, NAPROC
!/T JSEA = 1 + (ISEA-1)/NAPROC
!/T ETOT = 0.
!/T DO IK=1, NK
!/T E1I = 0.
!/T DO ITH=1, NTH
!/T E1I = E1I + A(ITH,IK,JSEA)
!/T END DO
!/T ETOT = ETOT + E1I * DSIP(IK) * SIG(IK) / CG(IK,ISEA)
!/T END DO
!/T IX = MAPSF(ISEA,1)
!/T IY = MAPSF(ISEA,2)
!/T HSIG (IX,IY) = 4. * SQRT ( ETOT * DTH )
!/T MAPOUT(IX,IY) = 1
!/T END DO
!
!/T IX0 = 1
!/T DO
!/T IXN = MIN ( NX , IX0+NXP-1 )
!/T CALL PRTBLK (NDST, NX, NY, NX, HSIG, MAPOUT, 0, 0., &
!/T IX0, IXN, 1, 1, NY, 1, 'Hs', 'm')
!/T IF ( IXN .EQ. NX ) EXIT
!/T IX0 = IX0 + NXP
!/T END DO
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3UINI : XGR = ',E10.3)
!
!/T1 9010 FORMAT (' TEST W3UINI : ISEA, U10C, XSTAR, ALPHA, FP, GAMMA')
!/T1 9011 FORMAT (' ',I6,F8.2,F10.1,2F6.3,F6.2)
!/
!/ End of W3UINI ----------------------------------------------------- /
!/
END SUBROUTINE W3UINI
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UBPT
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 12-Jun-2012 |
!/ +-----------------------------------+
!/
!/ 19-Oct-1998 : Final FORTRAN 77 ( version 1.18 )
!/ 20-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 07-Sep-2005 : Moving update to end of time step. ( version 3.08 )
!/ 17-Aug-2010 : Add initialization ABPI0-N(:,0). ( version 3.14.5 )
!/ 12-Jun-2012 : Add /RTD option or rotated grid option.
!/ (Jian-Guo Li) ( version 4.06 )
!/
! 1. Purpose :
!
! Update spectra at the active boundary points.
!
! 2. Method :
!
! Spectra are read and interpolated in space and time from the
! data read by W3IOBC.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
! STRACE, DSEC21
! Service routines.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - The data arrays contain sigma spectra to assure conservation
! when changing grids.
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T0 Test output of wave heights.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NSPEC, MAPWN, SIG2, DDEN
!/RTD !! Use rotation angle and action conversion sub. JGLi12Jun2012
!/RTD USE W3GDATMD, ONLY: NK, NTH, NSPEC, AnglD
!/RTD USE W3SERVMD, ONLY: w3acturn
USE W3ADATMD, ONLY: CG
USE W3ODATMD, ONLY: NBI, ABPI0, ABPIN, ISBPI, IPBPI, RDBPI, &
BBPI0, BBPIN
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: IBI, ISP, ISEA
!/S INTEGER, SAVE :: IENT = 0
!/T0 REAL :: HS1, HS2
!/RTD !! Declare a temporary spectr variable. JGLi12Jun2012
!/RTD REAL :: Spectr(NSPEC), AnglBP
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UBPT')
!
! 1. Process BBPI0 -------------------------------------------------- *
! 1.a First intialization
!
IF ( BBPI0(1,0) .EQ. -1. ) THEN
!
BBPI0(:,0) = 0.
BBPIN(:,0) = 0.
ABPI0(:,0) = 0.
ABPIN(:,0) = 0.
!
DO IBI=1, NBI
ISEA = ISBPI(IBI)
DO ISP=1, NSPEC
BBPI0(ISP,IBI) = CG(MAPWN(ISP),ISEA) / SIG2(ISP) * &
( RDBPI(IBI,1) * ABPI0(ISP,IPBPI(IBI,1)) &
+ RDBPI(IBI,2) * ABPI0(ISP,IPBPI(IBI,2)) &
+ RDBPI(IBI,3) * ABPI0(ISP,IPBPI(IBI,3)) &
+ RDBPI(IBI,4) * ABPI0(ISP,IPBPI(IBI,4)) )
END DO
END DO
!
! 1.b Shift BBPIN
!
ELSE
BBPI0 = BBPIN
END IF
!
! 2. Process BBPIN -------------------------------------------------- *
!
DO IBI=1, NBI
ISEA = ISBPI(IBI)
DO ISP=1, NSPEC
BBPIN(ISP,IBI) = CG(MAPWN(ISP),ISEA) / SIG2(ISP) * &
( RDBPI(IBI,1) * ABPIN(ISP,IPBPI(IBI,1)) &
+ RDBPI(IBI,2) * ABPIN(ISP,IPBPI(IBI,2)) &
+ RDBPI(IBI,3) * ABPIN(ISP,IPBPI(IBI,3)) &
+ RDBPI(IBI,4) * ABPIN(ISP,IPBPI(IBI,4)) )
END DO
!
!/RTD !! Rotate the spectra if model is on rotated grid. JGLi12Jun2012
!/RTD Spectr = BBPIN(:,IBI)
!/RTD AnglBP = AnglD(ISEA)
!/RTD CALL w3acturn( NTH, NK, AnglBP, Spectr )
!/RTD BBPIN(:,IBI) = Spectr
!/RTD
!
END DO
! 3. Wave height test output ---------------------------------------- *
!
!/T0 WRITE (NDST,9000)
!/T0 DO IBI=1, NBI
!/T0 HS1 = 0.
!/T0 HS2 = 0.
!/T0 DO ISP=1, NSPEC
!/T0 HS1 = HS1 + BBPI0(ISP,IBI) * DDEN(MAPWN(ISP)) / &
!/T0 CG(MAPWN(ISP),ISBPI(IBI))
!/T0 HS2 = HS2 + BBPIN(ISP,IBI) * DDEN(MAPWN(ISP)) / &
!/T0 CG(MAPWN(ISP),ISBPI(IBI))
!/T0 END DO
!/T0 HS1 = 4. * SQRT ( HS1 )
!/T0 HS2 = 4. * SQRT ( HS2 )
!/T0 WRITE (NDST,9001) IBI, ISBPI(IBI), HS1, HS2
!/T0 END DO
!
RETURN
!
! Formats
!
!/T0 9000 FORMAT ( ' TEST W3UBPT : WAVE HEIGHTS BBPI0/N (NO TAIL)')
!/T0 9001 FORMAT ( ' ',2I8,2X,2F8.2)
!/
!/ End of W3UBPT ----------------------------------------------------- /
!/
END SUBROUTINE W3UBPT
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UIC1( FLFRST )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | C. Sevigny |
!/ | FORTRAN 90 |
!/ | Last update : 27-Aug-2015 |
!/ +-----------------------------------+
!/
!/ 27-Aug-2015 : Creation ( version 5.10 )
!/
! 1. Purpose :
!
! Update ice thickness in the wave model.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! FLFRST L. I Spectra in 1-D or 2-D representation
! (points to same address).
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/DIST Id.
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NSEA, NSEA, MAPSF
USE W3WDATMD, ONLY: TIME, TIC1, ICEH
USE W3IDATMD, ONLY: TI1, ICEP1, FLIC1
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
LOGICAL, INTENT(IN) :: FLFRST
!/
!/ ------------------------------------------------------------------- /
!/ Local variables
!/
INTEGER :: IX, IY, ISEA
!/
!/
! 1. Preparations --------------------------------------------------- *
! 1.a Update times
!
!/T WRITE (NDST,9010) TIME, TIC1, TI1
TIC1(1) = TI1(1)
TIC1(2) = TI1(2)
! 2. Main loop over sea points -------------------------------------- *
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
ICEH(ISEA) = ICEP1(IX,IY)
END DO
!
RETURN
!/T 9010 FORMAT ( ' TEST W3UIC1 : TIME :',I9.8,I7.6/ &
!/T ' OLD TICE :',I9.8,I7.6/ &
!/T ' NEW TICE :',I9.8,I7.6)
!/
!/ End of W3UIC1 ----------------------------------------------------- /
!/
END SUBROUTINE W3UIC1
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UIC5( FLFRST )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | C. Sevigny & F. Ardhuin |
!/ | FORTRAN 90 |
!/ | Last update : 13-Jan-2016 |
!/ +-----------------------------------+
!/
!/ 27-Aug-2015 : Creation ( version 5.08 )
!/ 13-Jan-2016 : Changed initial value of ICEDMAX ( version 5.08 )
!/
! 1. Purpose :
!
! Update ice floe mean and max diameters in the wave model.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! FLFRST L. I Spectra in 1-D or 2-D representation
! (points to same address).
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/DIST Id.
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3IDATMD, ONLY: TI5, ICEP5
USE W3GDATMD, ONLY: NSEA, MAPSF
USE W3WDATMD, ONLY: TIME, TIC5, ICE, ICEH, ICEF, ICEDMAX
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
LOGICAL, INTENT(IN) :: FLFRST
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local variables
!/
INTEGER :: IX, IY, ISEA
LOGICAL :: FLFLOE
!/
!/
! 1. Preparations --------------------------------------------------- *
! 1.a Update times
!
!/T WRITE (NDST,9010) TIME, TIC5, TI5
TIC5(1) = TI5(1)
TIC5(2) = TI5(2)
! 2. Main loop over sea points -------------------------------------- *
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
FLFLOE = ICE(ISEA) .EQ. 0 .OR. ICEH(ISEA) .EQ. 0
IF ( FLFLOE) THEN
ICEF(ISEA) = 0.0
ICEDMAX(ISEA) = 1000.0
ELSE
ICEF(ISEA) = ICEP5(IX,IY)
ICEDMAX(ISEA) = ICEP5(IX,IY)
END IF
END DO
!
RETURN
!/T 9010 FORMAT ( ' TEST W3UIC5 : TIME :',I9.8,I7.6/ &
!/T ' OLD TICE :',I9.8,I7.6/ &
!/T ' NEW TICE :',I9.8,I7.6)
!/
!/
!/ End of W3UIC5 ----------------------------------------------------- /
!/
END SUBROUTINE W3UIC5
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UICE ( A, VA )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 28-Mar-2014 |
!/ +-----------------------------------+
!/
!/ 19-Oct-1998 : Final FORTRAN 77 ( version 1.18 )
!/ 20-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 11-Jan-2002 : Sub-grid ice. ( version 2.15 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 28-Jun-2005 : Adding MAPST2. ( version 3.07 )
!/ Taking out initilization.
!/ 11-May-2007 : Adding NTPROC/NAPROC separation. ( version 3.11 )
!/ 15-May-2010 : Adding second field for icebergs ( version 3.14 )
!/ 13-Mar-2012 : Add initialization of UST on re- ( version 4.07 )
!/ activation of grid point.
!/ 06-Jun-2012 : Porting bugfixes from 3.14 to 4.07 ( version 4.07 )
!/ 28-Mar-2014 : Adapting to ICx source terms ( version 4.18 )
!/
! 1. Purpose :
!
! Update ice map in the wave model.
!
! 2. Method :
!
! Points with an ice concentration larger than FICEN are removed
! from the sea map in the wave model. Such points are identified
! by negative numbers is the grid status map MAPSTA. For ice
! points spectra are set to zero. Points from wich ice disappears
! are initialized with a "small" JONSWAP spectrum, based on the
! frequency SIG(NK-1) and the local wind direction.
!
! In the case of icebergs, the iceberg attenuation coefficient is
! added to the subgrid obstruction map.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! (V)A R.A. I/O Spectra in 1-D or 2-D representation
! (points to same address).
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/DIST Id.
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF, MAPSTA, MAPST2, &
NTH, NK, NSPEC, SIG, TH, DTH, FICEN
USE W3WDATMD, ONLY: TIME, TICE, ICE, BERG, UST
!! USE W3ADATMD, ONLY: U10, U10D, CG
USE W3ADATMD, ONLY: CG
USE W3IDATMD, ONLY: TIN, ICEI, BERGI
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(INOUT) :: A(NTH,NK,0:NSEAL), VA(NSPEC,0:NSEAL)
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: IK, ITH, ISEA, JSEA, IX, IY, ISP
!/S INTEGER, SAVE :: IENT = 0
INTEGER :: MAPICE(NY,NX)
LOGICAL :: LOCAL
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UICE')
!
LOCAL = IAPROC .LE. NAPROC
!
!/T WRITE (NDST,9000) FICEN
!/T IF ( .NOT. LOCAL ) WRITE (NDST,9001)
!
! 1. Preparations --------------------------------------------------- *
! 1.a Update times
!
!/T WRITE (NDST,9010) TIME, TICE, TIN
TICE(1) = TIN(1)
TICE(2) = TIN(2)
!
! 1.b Process maps
!
!/IC0 MAPICE = MOD(MAPST2,2)
!/IC0 MAPST2 = MAPST2 - MAPICE
!
! 2. Main loop over sea points -------------------------------------- *
!
DO ISEA=1, NSEA
!
! 2.a Get grid counters
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
ICE(ISEA) = ICEI(IX,IY)
BERG(ISEA)= BERGI(IX,IY)
!
! 2.b Sea point to be de-activated..
!
!/IC0 IF ( ICEI(IX,IY).GE.FICEN .AND. MAPICE(IY,IX).EQ.0 ) THEN
!/IC0 MAPSTA(IY,IX) = - ABS(MAPSTA(IY,IX))
!/IC0 MAPICE(IY,IX) = 1
!/IC0 IF ( LOCAL .AND. MOD(ISEA-IAPROC,NAPROC) .EQ. 0 ) THEN
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'ICE (NEW)'
!/IC0 JSEA = 1 + (ISEA-1)/NAPROC
!/IC0 VA(:,JSEA) = 0.
!/IC0 ELSE
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'ICE (NEW X)'
!/IC0 END IF
!
!/IC0 ELSE IF ( ICEI(IX,IY).GE.FICEN ) THEN
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'ICE'
!/IC0 END IF
!
! 2.b Ice point to be re-activated.
!
!/IC0 IF ( ICEI(IX,IY).LT.FICEN .AND. MAPICE(IY,IX).EQ.1 ) THEN
!
!/IC0 MAPICE(IY,IX) = 0
!/IC0 UST(ISEA) = 0.05
!
!/IC0 IF ( MAPST2(IY,IX) .EQ. 0 ) THEN
!/IC0 MAPSTA(IY,IX) = ABS(MAPSTA(IY,IX))
!
!/IC0 IF ( LOCAL .AND. MOD(ISEA-IAPROC,NAPROC) .EQ. 0 ) THEN
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'SEA (NEW)'
!/IC0 JSEA = 1 + (ISEA-1)/NAPROC
!/IC0 VA(:,JSEA) = 0.
!
!/IC0 ELSE
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'SEA (NEW X)'
!/IC0 END IF
!
!/IC0 ELSE
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'DIS'
!/IC0 END IF
!
!/IC0 ELSE IF ( ICEI(IX,IY).LT.FICEN ) THEN
!/T WRITE (NDST,9021) ISEA, IX, IY, MAPSTA(IY,IX), &
!/T ICEI(IX,IY), 'SEA'
!
!/IC0 END IF
!
END DO
!
! 3. Update MAPST2 -------------------------------------------------- *
!
!/IC0 MAPST2 = MAPST2 + MAPICE
!
RETURN
!
! Formats
!
!/T 9000 FORMAT ( ' TEST W3UICE : FICEN :',F9.3)
!/T 9001 FORMAT ( ' TEST W3UICE : NO LOCAL SPECTRA')
!
!/T 9010 FORMAT ( ' TEST W3UICE : TIME :',I9.8,I7.6/ &
!/T ' OLD TICE :',I9.8,I7.6/ &
!/T ' NEW TICE :',I9.8,I7.6)
!
!/T 9020 FORMAT ( ' TEST W3UICE : ISEA, IX, IY, MAP, ICE, STATUS :')
!/T 9021 FORMAT ( ' ',I8,3I4,F6.2,2X,A)
!/
!/ End of W3UICE ----------------------------------------------------- /
!/
END SUBROUTINE W3UICE
!/ ------------------------------------------------------------------- /
SUBROUTINE W3ULEV ( A, VA )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 26-Sep-2012 |
!/ +-----------------------------------+
!/
!/ 15-Jan-1998 : Final FORTRAN 77 ( version 1.18 )
!/ 21-Jan-2000 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 30-Apr-2002 : Water level fixes. ( version 2.20 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 15-Jul-2005 : Adding drying out of points. ( version 3.07 )
!/ 11-May-2007 : Adding NTPROC/NAPROC separation. ( version 3.11 )
!/ 23-Aug-2011 : Bug fix for UG grids : new boundary ( version 4.04 )
!/ 13-Mar-2012 : Add initialization of UST on re- ( version 4.07 )
!/ activation of grid point.
!/ 06-Jun-2012 : Porting bugfixes from 3.14 to 4.07 ( version 4.07 )
!/ 26-Sep-2012 : Adding update from tidal analysis ( version 4.08 )
!/
! 1. Purpose :
!
! Update the water level.
!
! 2. Method :
!
! The wavenumber grid is modified without modyfying the spectrum
! (conservative linear interpolation to new grid).
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! (V)A R.A. I/O 2-D and 1-D represetation of the spectra.
! ----------------------------------------------------------------
!
! Local variables
! ----------------------------------------------------------------
! KDMAX Real Deep water cut-off for kd.
! WNO R.A. Old wavenumbers.
! CGO R.A. Old group velocities.
! OWN R.A. Old wavenumber band width.
! DWN R.A. New wavenumber band width.
! TA R.A. Auxiliary spectrum.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - The grid is updated only if KDmin > KDMAX.
! - The grid is updated for inactive points too.
! - The local wavenumber bandwidth is DSIGMA/CG.
! - The local spectrum is updated only if the grid is updated,
! the grid point is not disabled (MAPST2) and if the change of
! the lowest wavenumber exceeds RDKMIN times the band width.
! - No spectral initialization for newly wet points.
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Basic test output.
! !/T2 Output of minimum relative depth per grid point.
! !/T3 Spectra before and after
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF, MAPSTA, MAPST2, &
ZB, DMIN, NK, NTH, NSPEC, SIG, DSIP, &
MAPWN, MAPTH, FACHFA, GTYPE, UNGTYPE, W3SETREF
USE W3WDATMD, ONLY: TIME, TLEV, WLV, UST
USE W3ADATMD, ONLY: CG, WN, DW
USE W3IDATMD, ONLY: TLN, WLEV
USE W3SERVMD, ONLY: EXTCDE
USE W3DISPMD, ONLY: WAVNU1
USE W3TRIAMD, ONLY: SETUGIOBP
USE W3TIMEMD
!/TIDE USE W3GDATMD, ONLY: YGRD
!/TIDE USE W3IDATMD, ONLY: FLLEVTIDE, WLTIDE, NTIDE
!/TIDE USE W3TIDEMD
!/T3 USE W3ARRYMD, ONLY: PRT2DS
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(INOUT) :: A(NTH,NK,0:NSEAL), VA(NSPEC,0:NSEAL)
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: ISEA, JSEA, IX, IY, IK, I1, I2, &
ISPEC, IK0, ITH
!/S INTEGER, SAVE :: IENT = 0
INTEGER :: MAPDRY(NY,NX)
REAL :: DWO(NSEA), KDCHCK, WNO(0:NK+1), &
CGO(0:NK+1), DEPTH, &
RDK, RD1, RD2, TA(NTH,NK), &
OWN(NK), DWN(NK)
REAL :: KDMAX = 4., RDKMIN = 0.05
!/T3 REAL :: OUT(NK,NTH)
LOGICAL :: LOCAL
!
!/TIDE INTEGER :: J
!/TIDE INTEGER(KIND=4) :: TIDE_KD0, INT24, INTDYS ! "Gregorian day constant"
!/TIDE REAL :: WLEVTIDE, TIDE_ARG, WLEVTIDE2(1)
!/TIDE REAL(KIND=8) :: d1,h,TIDE_HOUR,HH,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau
!/TIDE REAL :: FX(44),UX(44),VX(44)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3ULEV')
!
LOCAL = IAPROC .LE. NAPROC
!
!/T WRITE (NDST,9000) KDMAX, RDKMIN
!
! 1. Preparations --------------------------------------------------- *
! 1.a Check NK
!
IF ( NK .LT. 2 ) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE (NDSE,1000)
CALL EXTCDE ( 1 )
END IF
!
! 1.b Update times
!
!/T WRITE (NDST,9010) TIME, TLEV
TLEV = TLN
!/T WRITE (NDST,9011) TLEV
!
! 1.c Extract dry point map, and residual MAPST2
!
MAPDRY = MOD(MAPST2/2,2)
MAPST2 = MAPST2 - 2*MAPDRY
!
! 1.d Update water levels and save old
!
!/TIDE IF (FLLEVTIDE) THEN
!/TIDE! WRITE(6,*) 'TIME:',TIME
!/TIDE TIDE_HOUR = TIME2HOURS(TIME)
!/TIDE!
!/TIDE!* THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
!/TIDE!* AT THE MID POINT OF THE ANALYSIS PERIOD.
!/TIDE d1=TIDE_HOUR/24.d0
!/TIDE TIDE_KD0= 2415020
!/TIDE d1=d1-dfloat(TIDE_kd0)-0.5d0
!/TIDE call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
!/TIDE INT24=24
!/TIDE INTDYS=int((TIDE_HOUR+0.00001)/INT24)
!/TIDE HH=TIDE_HOUR-dfloat(INTDYS*INT24)
!/TIDE TAU=HH/24.D0+H-S
!/TIDE END IF
!/TIDE!
!/TIDE! ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
!/TIDE! TING THE LUNAR TIME TAU.
!/TIDE!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
DWO(ISEA) = DW(ISEA)
!
!/TIDE IF (FLLEVTIDE) THEN
!/TIDE! VUF should be updated only if latitude changes significantly ...
!/TIDE CALL SETVUF_FAST(h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau,YGRD(IY,IX),FX,UX,VX)
!/TIDE WLEVTIDE = WLTIDE(IX,IY,1,1)
!/TIDE !Verification
!/TIDE ! IF (ISEA.EQ.1) THEN
!/TIDE
!/TIDE TIDE_AMPC(1:NTIDE,1)=WLTIDE(IX,IY,1:NTIDE,1)
!/TIDE TIDE_PHG(1:NTIDE,1)=WLTIDE(IX,IY,1:NTIDE,2)
!/TIDE!
!/TIDE ! WRITE(991,'(A,F20.2,13F8.3)') 'TEST ISEA 0:', &
!/TIDE ! d1,H,S,TAU,pp,s,p,enp,dh,dpp,ds,dp,dnp,YGRD(IY,IX)
!/TIDE J=1
!/TIDE ! WRITE(991,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',IX,J,TIME,TIDE_HOUR, &
!/TIDE ! FX(J),UX(J),VX(J),TIDE_INDEX2(J),TIDECON_ALLNAMES(TIDE_INDEX2(J))
!/TIDE DO J=2,TIDE_MF
!/TIDE TIDE_ARG=(VX(J)+UX(J))*twpi-WLTIDE(IX,IY,J,2)*DERA
!/TIDE WLEVTIDE =WLEVTIDE+FX(J)*WLTIDE(IX,IY,J,1)*COS(TIDE_ARG)
!/TIDE ! WRITE(991,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',IX,J,TIME,TIDE_HOUR, &
!/TIDE ! FX(J),UX(J),VX(J),TIDE_INDEX2(J),TIDECON_ALLNAMES(TIDE_INDEX2(J))
!/TIDE END DO
!/TIDE DO J=1,TIDE_MF
!/TIDE ! WRITE(991,'(A,4I9,F12.0,5F8.3)') 'TEST ISEA 2:',IX,J,TIME,TIDE_HOUR, &
!/TIDE ! FX(J),UX(J),VX(J),TIDE_AMPC(J,1),TIDE_PHG(J,1)
!/TIDE END DO
!/TIDE ! WRITE(991,'(A,3F7.3)') '#:',WLEV(IX,IY),WLEVTIDE,WLEV(IX,IY)-WLEVTIDE
!/TIDE ! CLOSE(991)
!/TIDE ! END IF
!/TIDE ! End of verification
!/TIDE WLV(ISEA) = WLEVTIDE
!/TIDE ELSE
!
WLV(ISEA) = WLEV(IX,IY)
!/TIDE ENDIF
DW (ISEA) = MAX ( 0. , WLV(ISEA)-ZB(ISEA) )
END DO
!
! 2. Loop over all sea points --------------------------------------- *
!
!/T2 WRITE (NDST,9020)
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
!
! 2.a Check if deep water
!
KDCHCK = WN(1,ISEA) * MIN( DWO(ISEA) , DW(ISEA) )
IF ( KDCHCK .LT. KDMAX ) THEN
!
! 2.b Update grid and save old grid
!
DEPTH = MAX ( DMIN, DW(ISEA) )
!
DO IK=0, NK+1
WNO(IK) = WN(IK,ISEA)
CGO(IK) = CG(IK,ISEA)
!
! Calculate wavenumbers and group velocities.
CALL WAVNU1(SIG(IK),DEPTH,WN(IK,ISEA),CG(IK,ISEA))
!
END DO
!
DO IK=1, NK
OWN(IK) = DSIP(IK) / CGO(IK)
DWN(IK) = DSIP(IK) / CG(IK,ISEA)
END DO
!
! 2.c Process dry points
!
IF ( WLV(ISEA)-ZB(ISEA) .LE.0. ) THEN
IF ( MAPDRY(IY,IX) .EQ. 0 ) THEN
IF ( LOCAL .AND. MOD(ISEA-IAPROC,NAPROC).EQ.0 ) THEN
JSEA = 1 + (ISEA-1)/NAPROC
VA(:,JSEA) = 0.
END IF
MAPDRY(IY,IX) = 1
MAPSTA(IY,IX) = -ABS(MAPSTA(IY,IX))
!/T2 WRITE (NDST,9021) ISEA, WLV(ISEA)-ZB(ISEA), &
!/T2 0., 0., ' (NEW DRY)'
!/T2 ELSE
!/T2 WRITE (NDST,9021) ISEA, WLV(ISEA)-ZB(ISEA), &
!/T2 0., 0., ' (DRY)'
ENDIF
CYCLE
END IF
!
! 2.d Process new wet point
!
IF (WLV(ISEA)-ZB(ISEA).GT.0. .AND. MAPDRY(IY,IX).EQ.1) THEN
MAPDRY(IY,IX) = 0
!
! Resets the spectrum to zero
!
IF ( LOCAL .AND. MOD(ISEA-IAPROC,NAPROC).EQ.0 ) THEN
JSEA = 1 + (ISEA-1)/NAPROC
VA(:,JSEA) = 0.
END IF
!
UST(ISEA) = 0.05
IF ( MAPST2(IY,IX) .EQ. 0 ) THEN
MAPSTA(IY,IX) = ABS(MAPSTA(IY,IX))
!/T2 WRITE (NDST,9021) ISEA, WLV(ISEA)-ZB(ISEA), &
!/T2 0., 0., ' (NEW WET)'
!/T2 ELSE
!/T2 WRITE (NDST,9021) ISEA, WLV(ISEA)-ZB(ISEA), &
!/T2 0., 0., ' (NEW WET INACTIVE)'
END IF
CYCLE
END IF
!
! 2.e Check if ice on grid point, or if grid changes negligible
!
RDK = ABS(WNO(1)-WN(1,ISEA)) / DWN(1)
!
!/T2 IF ( MAPSTA(IY,IX) .LT. 0 ) THEN
!/T2 WRITE (NDST,9021) &
!/T2 ISEA, DW(ISEA), KDCHCK, RDK, ' (INACTIVE)'
!/T2 ELSE IF ( RDK .LT. RDKMIN ) THEN
!/T2 WRITE (NDST,9021) &
!/T2 ISEA, DW(ISEA), KDCHCK, RDK, ' (NEGL)'
!/T2 ELSE
!/T2 WRITE (NDST,9021) &
!/T2 ISEA, DW(ISEA), KDCHCK, RDK, ' '
!/T2 END IF
!
IF ( RDK.LT.RDKMIN .OR. MAPSTA(IY,IX).LT.0 ) CYCLE
IF ( MOD(ISEA-IAPROC,NAPROC) .NE. 0 ) CYCLE
JSEA = 1 + (ISEA-1)/NAPROC
!
IF ( .NOT. LOCAL ) CYCLE
!
! 2.d Save discrete actions and clean spectrum
!
DO IK=1, NK
DO ITH=1, NTH
!/T3 OUT(IK,ITH) = A(ITH,IK,JSEA) * SIG(IK) / CGO(IK)
TA(ITH,IK) = A(ITH,IK,JSEA) * OWN(IK)
END DO
END DO
!
VA(:,JSEA) = 0.
!
!/T3 CALL PRT2DS ( NDST, NK, NK, NTH, OUT, SIG(1), ' ', &
!/T3 TPI, 0., 1.E-5, 'F(f,th)', 'm2s', 'Before' )
!
! 2.e Redistribute discrete action density
!
IF ( WNO(1) .LT. WN(1,ISEA) ) THEN
IK0 = 1
I1 = 0
I2 = 1
220 CONTINUE
IK0 = IK0 + 1
IF ( IK0 .GT. NK+1 ) GOTO 251
IF ( WNO(IK0) .GE. WN(1,ISEA) ) THEN
IK0 = IK0 - 1
ELSE
GOTO 220
END IF
ELSE
IK0 = 1
I1 = 1
I2 = 2
END IF
!
DO 250, IK=IK0, NK
!
230 CONTINUE
IF ( WNO(IK) .GT. WN(I2,ISEA) ) THEN
I1 = I1 + 1
IF ( I1 .GT. NK ) GOTO 250
I2 = I1 + 1
GOTO 230
END IF
!
IF ( I1 .EQ. 0 ) THEN
RD1 = ( WN(1,ISEA) - WNO(IK) ) / DWN(1)
RD2 = 1. - RD1
ELSE
RD1 = ( WN(I2,ISEA) - WNO(IK) ) / &
( WN(I2,ISEA) - WN(I1,ISEA) )
RD2 = 1. - RD1
END IF
!
IF ( I1 .GE. 1 ) THEN
DO ITH=1, NTH
A(ITH,I1,JSEA) = A(ITH,I1,JSEA) + RD1*TA(ITH,IK)
END DO
END IF
!
IF ( I2 .LE. NK ) THEN
DO ITH=1, NTH
A(ITH,I2,JSEA) = A(ITH,I2,JSEA) + RD2*TA(ITH,IK)
END DO
END IF
!
250 CONTINUE
251 CONTINUE
!
! 2.f Convert discrete action densities to spectrum
!
DO ISPEC=1, NSPEC
VA(ISPEC,JSEA) = VA(ISPEC,JSEA) / DWN(MAPWN(ISPEC))
END DO
!
! 2.f Add tail if necessary
!
IF ( I2.LE.NK .AND. RD2.LE.0.95 ) THEN
DO IK=MAX(I2,2), NK
DO ITH=1, NTH
A(ITH,IK,JSEA) = FACHFA * A(ITH,IK-1,JSEA)
END DO
END DO
END IF
!
!/T3 DO ISPEC=1, NSPEC
!/T3 IK = MAPWN(ISPEC)
!/T3 ITH = MAPTH(ISPEC)
!/T3 OUT(IK,ITH) = A(ITH,IK,JSEA) * SIG(IK) / CG(IK,ISEA)
!/T3 END DO
!
!/T3 CALL PRT2DS ( NDST, NK, NK, NTH, OUT, SIG(1), ' ', &
!/T3 TPI, 0., 1.E-5, 'F(f,th)', 'm2s', 'After' )
!
!/T2 ELSE
!/T2 WRITE (NDST,9021) ISEA, KDCHCK, ' (DEEP)'
END IF
!
END DO
!
! 3. Reconstruct new MAPST2 ----------------------------------------- *
!
MAPST2 = MAPST2 + 2*MAPDRY
!
! 4. Re-generates the boundary data ---------------------------------- *
!
IF (GTYPE.EQ.UNGTYPE) THEN
CALL SETUGIOBP
!/REF1 ELSE
!/REF1 CALL W3SETREF
ENDIF
!
RETURN
!
! Formats
!
1000 FORMAT (/' *** ERROR W3ULEV *** '/ &
' THIS ROUTINE REQUIRES NK > 1 '/)
!
!/T 9000 FORMAT ( ' TEST W3ULEV : KDMAX :',F6.1/ &
!/T ' RDKMIN :',F8.3)
!
!/T 9010 FORMAT ( ' TEST W3ULEV : TIME :',I9.8,I7.6/ &
!/T ' OLD TLEV :',I9.8,I7.6)
!/T 9011 FORMAT ( ' NEW TLEV :',I9.8,I7.6)
!
!/T2 9020 FORMAT ( ' TEST W3ULEV : LOOP OVER ALL POINTS:', &
!/T2 ' ISEA, DW, KDMIN, RDK : ')
!/T2 9021 FORMAT ( ' ',I6,F8.2,F6.2,F7.3,A)
!/
!/ End of W3ULEV ----------------------------------------------------- /
!/
END SUBROUTINE W3ULEV
!/ ------------------------------------------------------------------- /
SUBROUTINE W3UTRN ( TRNX, TRNY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 30-Oct-2009 |
!/ +-----------------------------------+
!/
!/ 02-Apr-2001 : Origination. ( version 2.10 )
!/ 11-Jan-2002 : Sub-grid ice. ( version 2.15 )
!/ 30-Apr-2002 : Change to ICE on storage grid. ( version 2.20 )
!/ 15-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 11-Jan-2007 : Clean-up for boundary points. ( version 3.10 )
!/ 30-Oct-2009 : Implement run-time grid selection. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/
! 1. Purpose :
!
! Update cell boundary transparencies for general use in propagation
! routines.
!
! 2. Method :
!
! Two arrays are generated with the size (NY*NX,-1:1). The value
! at (IXY,-1) indicates the transparency to be used if the lower
! or left boundary is an inflow boundary. (IXY,1) is used if the
! upper or right boundary is an inflow boundary. (IXY,0) is used
! for all other cases (by definition full transparency).
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! TRNX/Y R.A. I Transparencies from model defintion file.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Basic test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSTA, MAPSF, &
TRFLAG, FICE0, FICEN, FICEL, &
RLGTYPE, CLGTYPE, GTYPE, FLAGLL, &
HPFAC, HQFAC, FFACBERG
USE W3WDATMD, ONLY: ICE, BERG
USE W3ADATMD, ONLY: ATRNX, ATRNY
!
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: TRNX(NY*NX), TRNY(NY*NX)
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: ISEA, IX, IY, IXY, IXN, IXP, IYN, IYP
!/S INTEGER, SAVE :: IENT = 0
!/T INTEGER :: ILEV, NLEV
REAL :: TRIX(NY*NX), TRIY(NY*NX), DX, DY, &
LICE0, LICEN
!/T REAL :: LEVS(0:10)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3UTRN')
!/T WRITE (NDST,9000) TRFLAG
!
! 1. Preparations --------------------------------------------------- *
!
ATRNX = 1.
ATRNY = 1.
!/T WRITE (NDST,9001) 'INITIALIZING ATRNX/Y'
!
! 2. Filling arrays from TRNX/Y for obstructions -------------------- *
! 2.a TRFLAG = 0, no action needed
IF ( TRFLAG .EQ. 0 ) THEN
!/T WRITE (NDST,9001) 'NO FURTHER ACTION REQUIRED'
RETURN
!
! 2.b TRFLAG = 1,3: TRNX/Y defined at boundaries
!
ELSE IF ( TRFLAG.EQ.1 .OR. TRFLAG.EQ.3 .OR. TRFLAG.EQ.5 ) THEN
!/T WRITE (NDST,9001) 'DATA APPLIED AT CELL BOUNDARIES'
!/T LEVS = 0.
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
IF ( IX .EQ. 1 ) THEN
ATRNX(IXY,-1) = TRNX(IY+(NX-1)*NY)
ATRNX(IXY, 1) = TRNX(IXY)
ELSE IF ( IX .EQ. NX ) THEN
ATRNX(IXY,-1) = TRNX(IXY-NY)
ATRNX(IXY, 1) = TRNX(IY)
ELSE
ATRNX(IXY,-1) = TRNX(IXY-NY)
ATRNX(IXY, 1) = TRNX(IXY)
END IF
ATRNY(IXY,-1) = TRNY(IXY-1)
ATRNY(IXY, 1) = TRNY(IXY)
!
!/T ILEV = NINT(10.*MIN(TRNX(IXY),TRNY(IXY)))
!/T LEVS(ILEV) = LEVS(ILEV) + 1.
!
END DO
!
! 2.c TRFLAG = 2,4: TRNX/Y defined at cell centers
!
ELSE
!/T WRITE (NDST,9001) 'DATA APPLIED AT CELL CENTERS'
!/T LEVS = 0.
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
!
IF ( IX .EQ. 1 ) THEN
IXN = IY + (NX-1)*NY
IXP = IY + IX *NY
ELSE IF ( IX .EQ. NX ) THEN
IXN = IY + (IX-2)*NY
IXP = IY
ELSE
IXN = IY + (IX-2)*NY
IXP = IY + IX *NY
END IF
!
IF ( IY .EQ. 1 ) THEN
IYN = IXY
IYP = IXY + 1
ELSE IF ( IY .EQ. NY ) THEN
IYN = IXY - 1
IYP = IXY
ELSE
IYN = IXY - 1
IYP = IXY + 1
END IF
!
! factors 0.5 in first term and 2. in second term cancel
!
ATRNX(IXY,-1) = (1.+TRNX(IXY)) * TRNX(IXN)/(1.+TRNX(IXN))
ATRNX(IXY, 1) = (1.+TRNX(IXY)) * TRNX(IXP)/(1.+TRNX(IXP))
ATRNY(IXY,-1) = (1.+TRNY(IXY)) * TRNY(IYN)/(1.+TRNY(IYN))
ATRNY(IXY, 1) = (1.+TRNY(IXY)) * TRNY(IYP)/(1.+TRNY(IYP))
!
IF ( MAPSTA(IY,IX) .EQ. 2 ) THEN
IF ( IX .EQ. 1 ) THEN
ATRNX(IXY,-1) = 1.
ELSE IF ( MAPSTA( IY ,IX-1) .LE. 0 ) THEN
ATRNX(IXY,-1) = 1.
END IF
IF ( IX .EQ. NX ) THEN
ATRNX(IXY, 1) = 1.
ELSE IF ( MAPSTA( IY ,IX+1) .LE. 0 ) THEN
ATRNX(IXY, 1) = 1.
END IF
IF ( IY .EQ. 1 ) THEN
ATRNY(IXY,-1) = 1.
ELSE IF ( MAPSTA(IY-1, IX ) .LE. 0 ) THEN
ATRNY(IXY,-1) = 1.
END IF
IF ( IY .EQ. NY ) THEN
ATRNY(IXY, 1) = 1.
ELSE IF ( MAPSTA(IY+1, IX ) .LE. 0 ) THEN
ATRNY(IXY, 1) = 1.
END IF
END IF
!
!/T ILEV = NINT(10.*MIN(TRNX(IXY),TRNY(IXY)))
!/T LEVS(ILEV) = LEVS(ILEV) + 1.
!
END DO
END IF
!
!/T WRITE(NDST,9010) 'ISLANDS'
!/T NLEV = 0
!/T DO ILEV=0, 10
!/T WRITE (NDST,9011) ILEV, LEVS(ILEV)/REAL(NSEA)
!/T NLEV = NLEV + NINT(LEVS(ILEV))
!/T END DO
!
! 3. Adding ice to obstructions ------------------------------------- *
! 3.a TRFLAG < 3, no action needed
!
IF ( TRFLAG.LT.3 .OR. FICEN-FICE0.LT.1.E-6 ) THEN
!/T WRITE (NDST,9001) 'NO ICE ACTION REQUIRED'
RETURN
!
! 3.b TRFLAG = 3,4: Calculate ice transparencies
!
ELSE
!/T WRITE (NDST,9001) 'CALCULATE ICE TRANSPARENCIES'
!/T LEVS = 0.
TRIX = 1.
TRIY = 1.
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
!
DX = HPFAC(IY,IX)
DY = HQFAC(IY,IX)
IF ( FLAGLL ) THEN
DX = DX * RADIUS * DERA
DY = DY * RADIUS * DERA
END IF
!
!/IC0 IF (ICE(ISEA).GT.0) THEN
!/IC0 IF (FICEL.GT.0.) THEN
!/IC0 TRIX(IXY) = EXP(-ICE(ISEA)*DX/FICEL)
!/IC0 TRIY(IXY) = EXP(-ICE(ISEA)*DY/FICEL)
!/IC0 ELSE
! Otherwise: original Tolman expression (Tolman 2003)
!/IC0 LICE0 = FICE0*DX
!/IC0 LICEN = FICEN*DX
!/IC0 TRIX(IXY) = ( LICEN - ICE(ISEA)*DX ) / ( LICEN - LICE0 )
! begin temporary notes
! TRIX = ( LICEN - ICE(ISEA)*DX ) / ( LICEN - LICE0 )
! thus, it is TRIX= ( (FICEN*DX) - ICE(ISEA)*DX ) / ( (FICEN*DX) - (FICE0*DX) )
! thus, it is TRIX= ( FICEN - ICE(ISEA) ) / ( FICEN - FICE0 )
! in other words, the variables DX DY are not used
! and the variables LICE0 LICEN are not necessary.
! end temporary notes
!/IC0 LICE0 = FICE0*DY
!/IC0 LICEN = FICEN*DY
!/IC0 TRIY(IXY) = ( LICEN - ICE(ISEA)*DY ) / ( LICEN - LICE0 )
!/IC0 END IF
!
!/IC0 TRIX(IXY) = MAX ( 0. , MIN ( 1. , TRIX(IXY) ) )
!/IC0 TRIY(IXY) = MAX ( 0. , MIN ( 1. , TRIY(IXY) ) )
!/IC0 END IF
!
! Adding iceberg attenuation
!
IF (BERG(ISEA).GT.0) THEN
TRIX(IXY) = TRIX(IXY)*EXP(-BERG(ISEA)*FFACBERG *DX*0.0001)
TRIY(IXY) = TRIY(IXY)*EXP(-BERG(ISEA)*FFACBERG *DY*0.0001)
END IF
!
!/T ILEV = NINT(10.*MIN(TRIX(IXY),TRIY(IXY)))
!/T LEVS(ILEV) = LEVS(ILEV) + 1.
!
END DO
!
!/T WRITE(NDST,9010) 'ICE'
!/T NLEV = 0
!/T DO ILEV=0, 10
!/T WRITE (NDST,9011) ILEV, LEVS(ILEV)/REAL(NSEA)
!/T NLEV = NLEV + NINT(LEVS(ILEV))
!/T END DO
!
! 3.c Combine transparencies, ice always defined at cell center !
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
!
IF ( IX .EQ. 1 ) THEN
IXN = IY + (NX-1)*NY
IXP = IY + IX *NY
ELSE IF ( IX .EQ. NX ) THEN
IXN = IY + (IX-2)*NY
IXP = IY
ELSE
IXN = IY + (IX-2)*NY
IXP = IY + IX *NY
END IF
!
IF ( IY .EQ. 1 ) THEN
IYN = IXY
IYP = IXY + 1
ELSE IF ( IY .EQ. NY ) THEN
IYN = IXY - 1
IYP = IXY
ELSE
IYN = IXY - 1
IYP = IXY + 1
END IF
!
ATRNX(IXY,-1) = ATRNX(IXY,-1) &
* (1.+TRIX(IXY)) * TRIX(IXN)/(1.+TRIX(IXN))
ATRNX(IXY, 1) = ATRNX(IXY, 1) &
* (1.+TRIX(IXY)) * TRIX(IXP)/(1.+TRIX(IXP))
ATRNY(IXY,-1) = ATRNY(IXY,-1) &
* (1.+TRIY(IXY)) * TRIY(IYN)/(1.+TRIY(IYN))
ATRNY(IXY, 1) = ATRNY(IXY, 1) &
* (1.+TRIY(IXY)) * TRIY(IYP)/(1.+TRIY(IYP))
!
END DO
!
END IF
!
RETURN
!
! Formats
!
!/T 9000 FORMAT ( ' TEST W3UTRN : TRFLAG = ',I3)
!/T 9001 FORMAT ( ' TEST W3UTRN : ',A)
!/T 9010 FORMAT ( ' TEST W3UTRN : OBSTRICTION LEVELS FOR ',A,' :')
!/T 9011 FORMAT ( ' ',I4,F8.5)
!/
!/ End of W3UTRN ----------------------------------------------------- /
!/
END SUBROUTINE W3UTRN
!/ ------------------------------------------------------------------- /
SUBROUTINE W3DZXY( ZZ, ZUNIT, DZZDX, DZZDY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | W. E. Rogers, NRL |
!/ | FORTRAN 90 |
!/ | Last update : 06-Dec-2010 |
!/ +-----------------------------------+
!/
!/ 30-Oct-2009 : Origination. ( version 3.14 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/
! 1. Purpose :
!
! Calculate derivatives of a field.
!
! 2. Method :
!
! Derivatives are calculated in m/m from the longitude/latitude
! grid, central in space for iternal points, one-sided for
! coastal points.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ZZ R.A. I Field to calculate derivatives of.
! ZUNIT R.A. I Units of ZZ (used for test output).
! DZZDX R.A. O Derivative in X-direction (W-E).
! DZZDY R.A. O Derivative in Y-direction (S-N).
! IXP: IX plus 1 (with branch cut incorporated)
! IYP, IXM, IYM: ditto
! IXPS: value to use for IXP if IXPS is not masked.
! (use IX if masked)
! IYPS, IXMS, IYMS : ditto
! IXTRPL : in case of needing IY+1 for IY=NY, IX needs to be
! modified (tripole grid only)
! IXTRPLS : value to use for IXTRPL if IXTRPLS is not masked
! (use IX if masked)
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! This routine replaces the functionality of W3DDXY and W3DCXY.
!
! Output arrays are always initialized to zero.
!
! 8. Structure :
!
! ----------------------------------------
! 1. Preparations
! a Initialize arrays
! b Set constants
! 2. Derivatives in X-direction (W-E).
! 3. Derivatives in Y-direction (S-N).
! ----------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSTA, MAPFS, MAPFS, &
DPDX, DPDY, DQDX, DQDY, FLAGLL, ICLOSE, &
ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPERR, NAPROC
USE W3SERVMD, ONLY: EXTCDE
!/T USE W3ARRYMD, ONLY : PRTBLK
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
REAL, INTENT(IN) :: ZZ(NSEA)
CHARACTER, INTENT(IN) :: ZUNIT*(*)
REAL, INTENT(OUT) :: DZZDX(NY,NX), DZZDY(NY,NX)
INTEGER :: ISEA, IX, IY, IXP, IXM, IYP, IYM
!/T INTEGER :: ISX, ISY, MAPOUT(NX,NY)
!/S INTEGER, SAVE :: IENT = 0
!/T INTEGER, SAVE :: NXS = 49
REAL :: DFAC , STX, STY
INTEGER :: IXPS,IYPS,IXMS,IYMS,IXTRPL,IXTRPLS
INTEGER :: IXSTART,IXEND
!/T REAL :: XOUT(NX,NY)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3DZXY')
!
! 1. Preparations --------------------------------------------------- *
! 1.a Initialize arrays
!
DZZDX = 0.
DZZDY = 0.
!
! 1.b Set constants
!
IF ( FLAGLL ) THEN
DFAC = 1. / ( DERA * RADIUS )
ELSE
DFAC = 1.
END IF
!
! 2. Derivatives in X-direction (W-E) and Y-direction (S-N) ----- *
!
! 2a. All points previously done in 2a,2b,2c of v4.18 done in 2a now:
IF ( ICLOSE.EQ.ICLOSE_NONE ) THEN
IXSTART=2
IXEND=NX-1
ELSE
IXSTART=1
IXEND=NX
ENDIF
DO IY=2, NY-1
DO IX=IXSTART,IXEND
IF ( MAPSTA(IY,IX) .NE. 0 ) THEN
STX = 0.5
IF (IX.EQ.NX)THEN
IXPS=1
ELSE
IXPS=IX+1
ENDIF
IF (MAPSTA(IY,IXPS).EQ.0) THEN
IXP = IX
STX = 1.0
ELSE
IXP = IXPS
END IF
IF(IX.EQ.1)THEN
IXMS=NX
ELSE
IXMS=IX-1
ENDIF
IF (MAPSTA(IY,IXMS).EQ.0) THEN
IXM = IX
STX = 1.0
ELSE
IXM = IXMS
END IF
STY = 0.5
IYPS=IY+1
IF (MAPSTA(IYPS,IX).EQ.0) THEN
IYP = IY
STY = 1.0
ELSE
IYP = IYPS
END IF
IYMS=IY-1
IF (MAPSTA(IYMS,IX).EQ.0) THEN
IYM = IY
STY = 1.0
ELSE
IYM = IYMS
END IF
DZZDX(IY,IX) = (ZZ(MAPFS(IY ,IXP))-ZZ(MAPFS(IY ,IXM))) * STX * DPDX(IY,IX) &
+ (ZZ(MAPFS(IYP,IX ))-ZZ(MAPFS(IYM,IX ))) * STY * DQDX(IY,IX)
DZZDY(IY,IX) = (ZZ(MAPFS(IY ,IXP))-ZZ(MAPFS(IY ,IXM))) * STX * DPDY(IY,IX) &
+ (ZZ(MAPFS(IYP,IX ))-ZZ(MAPFS(IYM,IX ))) * STY * DQDY(IY,IX)
DZZDX(IY,IX) = DZZDX(IY,IX) * DFAC
DZZDY(IY,IX) = DZZDY(IY,IX) * DFAC
END IF
END DO
END DO
! 2b. column IY=NY for tripole case
! This is more complex, since for these two points: (IYP,IX) (IYM,IX),
! not only is the first index different (IYP.NE.IYM), but also the
! second index is different (IX.NE.IX)!
IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
IY=NY
DO IX=1, NX
IF ( MAPSTA(IY,IX) .NE. 0 ) THEN
STX = 0.5
IF (IX.EQ.NX)THEN
IXPS=1
ELSE
IXPS=IX+1
ENDIF
IF (MAPSTA(IY,IXPS).EQ.0) THEN
IXP = IX
STX = 1.0
ELSE
IXP = IXPS
END IF
IF(IX.EQ.1)THEN
IXMS=NX
ELSE
IXMS=IX-1
ENDIF
IF (MAPSTA(IY,IXMS).EQ.0) THEN
IXM = IX
STX = 1.0
ELSE
IXM = IXMS
END IF
STY = 0.5
!..............next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
!..............i.e. target point is MAPFS(IY,(NX-IX+1))
IXTRPLS=NX-IX+1
IF (MAPSTA(IY,IXTRPLS).EQ.0) THEN
IXTRPL = IX
STY = 1.0
ELSE
IXTRPL=IXTRPLS
END IF
IYMS=IY-1
IF (MAPSTA(IYMS,IX).EQ.0) THEN
IYM = IY
STY = 1.0
ELSE
IYM = IYMS
END IF
! tripole grid: (IYP,IX) is replaced with (IY,IXTRPL)
DZZDX(IY,IX) = (ZZ(MAPFS(IY ,IXP))-ZZ(MAPFS(IY ,IXM))) * STX * DPDX(IY,IX) &
+ (ZZ(MAPFS(IY,IXTRPL))-ZZ(MAPFS(IYM,IX ))) * STY * DQDX(IY,IX)
DZZDY(IY,IX) = (ZZ(MAPFS(IY ,IXP))-ZZ(MAPFS(IY ,IXM))) * STX * DPDY(IY,IX) &
+ (ZZ(MAPFS(IY,IXTRPL))-ZZ(MAPFS(IYM,IX ))) * STY * DQDY(IY,IX)
DZZDX(IY,IX) = DZZDX(IY,IX) * DFAC
DZZDY(IY,IX) = DZZDY(IY,IX) * DFAC
END IF
END DO
END IF ! IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
!
! 3. Test output of fields ------------------------------------------ *
!
!/T WRITE (NDST,9010)
!/T ISX = 1 + NX/NXS
!/T ISY = 1 + NY/NXS
!/T DO IY=1, NY
!/T DO IX=1, NX
!/T MAPOUT(IX,IY) = MAPSTA(IY,IX)
!/T IF ( MAPFS(IY,IX) .NE. 0 ) &
!/T XOUT(IX,IY) = ZZ(MAPFS(IY,IX))
!/T END DO
!/T END DO
!/T CALL PRTBLK (NDST, NX, NY, NX, XOUT, MAPOUT, 0, 0., &
!/T 1, NX, ISX, 1, NY, ISY, 'ZZ', ZUNIT)
!/T DO IY=1, NY
!/T DO IX=1, NX
!/T XOUT(IX,IY) = DZZDX(IY,IX)
!/T END DO
!/T END DO
!/T CALL PRTBLK (NDST, NX, NY, NX, XOUT, MAPOUT, 0, 0., &
!/T 1, NX, ISX, 1, NY, ISY, 'DZZDX',TRIM(ZUNIT)//'/m')
!/T DO IY=1, NY
!/T DO IX=1, NX
!/T XOUT(IX,IY) = DZZDY(IY,IX)
!/T END DO
!/T END DO
!/T CALL PRTBLK (NDST, NX, NY, NX, XOUT, MAPOUT, 0, 0., &
!/T 1, NX, ISX, 1, NY, ISY, 'DZZDY',TRIM(ZUNIT)//'/m')
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3DZXY : DX0I, DY0I : ',2E12.5)
!/T 9010 FORMAT (' TEST W3DZXY : FIELDS ')
!/
!/ End of W3DZXY ----------------------------------------------------- /
!/
END SUBROUTINE W3DZXY
!/ ------------------------------------------------------------------- /
!/ End of module W3UPDTMD -------------------------------------------- /
!/
END MODULE W3UPDTMD