Definition of depths in WW3
The depths read by ww3_grid to generate the configuration are by default defined positive up (negative under water). There are two depths threshols: ZLIM is the value above which nodes will be considered as "forever dry land", and DMIN is the minimum water depth in the computation.
For example, if you have a bathymetry file with values -8 -2 2 and that you set ZLIM=0.5 and DMIN=2 then the model will compute the waves at the first two nodes (-8 and -2) but not at the one with Z = 2, because this is > 0.5 . If during the computation you add a water level (e.g. tidal elevation) of LEV=+3 m, the computation will still not be done at the last point because what matters is the values in the depth file compared to ZLIM.
If you add LEV = -1, then the water depth at the second point becomes DEPTH = -1 - (-2 ) = 1 m, but the model will treat it as if it was 2 m: this is because of DMIN = 2. The only place where the depth is not set to DMIN is the Miche limiter (which was buggy but just got fixed in version 4.12) and which uses the true depth, not DMIN.
If you add LEV = -3 (a lower water level) this is point becomes dry since -3 - (-2 ) = -1 < 0.
Wave direction definition
This one is really mind-boggling: the directions in the code are given by the variable THETA, and they represent the directions TO which the waves are propagating, with a trigonometric convention (0 = to the east). These are in radians and take values from 0 to 2*pi*(NTH-1)/NTH . OK. Now, on output, for example in the spectral output, the angles are still ordered in the same way (first direction is to East ...) but the values given in the angle vector have been converted to nautical convention ... BUT ... in radians. This is why these numbers start at pi/2, decrease to 0, then jump to 2*pi*(NTH-1)/NTH and decrease again to pi/2+2*p/NTH ... easy, uh?
Here is what it gives for 24 directions: this is in the header of the *spc files generated by ww3_outp :
0.157E+01 0.131E+01 0.105E+01 0.785E+00 0.524E+00 0.262E+00 0.000E+00 0.602E+01 0.576E+01 0.550E+01 0.524E+01 0.497E+01 0.471E+01 0.445E+01 0.419E+01 0.393E+01 0.367E+01 0.340E+01 0.314E+01 0.288E+01 0.262E+01 0.236E+01 0.209E+01 0.183E+01
What it means is (in degrees, nautical convention)
90 (= to east) 75 60 45 30 15 0 (to north) -15 -30 -45 ....
WWATCH is built around calls to W3WAVE which performs an integration of the model over a given time frame. Once this time integration is performed, W3WAVE finishes and WW3_SHEL determines the next time frame and calls W3WAVE again... this is repeated until the final time TIMEN is reached.
Program WW3_SHEL calls subroutine W3WAVE several times until TIMEN is reached. In a coupled system WW3_SHEL is itself a subroutine.
- The current time is given by the variable TIME
- Time in WAVEWATCH is defined with an INTEGER(2) array. For example TIME(1)=20120121 and TIME(2)=133421 means that TIME corresponds to Jan 21st 2012 at 1 PM, 34 min and 21 sec. Note that with this system, WWATCH cannot tell the difference between times closer than 1 second.
- specific routines for comparing different times or incrementing the clock are contained in the module W3TIMEMD.
- The time of the run go from TIME0 to TIMEN. Be careful that after being read from ww3_shel.inp, TIME0 is redefined and then stands for the end of the next integration step.
- Input times are defined by TFN(1:2,J), where the first two indices give the date or hour, and J is the input type.
- Input fields (currents and winds) are updated by a call to W3UCUR or W3UWND inside of W3WAVE. Be careful that winds and currents can be interpolated in time, which requires the knowledge of the current or wind at a future time level. This may not fly with a coupler... instead the coupler should probably provide directly the CX and CY fields (U and V current components) which correspond the the fields after time interpolation by W3UCUR.
- The TFN array is stored in the W3IDATMD module, and its values are updated by W3FLD??
- Output (and coupling) time steps are defined by the TONEXT array in the loop J = 1, NOTYPE of ww3_shel. This array is held in the module W3ODATMD and is defined by the routine W3INIT which is in the module WMINITMD .
- DSEC21 ( TTT , TTIME ) gives back the difference (in seconds) between the two times TTT and TTIME. This is used in tests to define the next input time step, leading to the definition of TIME0 which is the next time until which the model is integrated: CALL W3WAVE ( 1, TIME0 )
The top-level time integration loop starts at 700 CONTINUE in the program WW3_SHEL and finishes at GOTO 700. Inside this loop you will see that WWATCH first defines the integration time frame all the way to the finish time TIMEN. This is done with
TTIME = TIMEN
However, after this bold attempt, WWATCH then checks on all the possible input times in between and accordingly resets TTIME to the next input time with this statement
TTIME = TTT
But what about output ??? It looks like the output times are only used inside the W3WAVE subroutine. Indeed in W3WAVE you will find this :
DO J=1, NOTYPE IF ( FLOUT(J) ) THEN ! ! 4.d Perform output ! TOUT(:) = TONEXT(:,J) DTTST = DSEC21 ( TIME, TOUT )
Which looks like several output (and thus coupling) steps can be done without going out of W3WAVE back to WW3_SHEL. Hence for coupling WWATCH to another model, if we want to have the coupler calls outside of W3WAVE we have to make sure that the TTIME in WW3_SHEL is incremented at times when an export to the coupler is needed. Assuming that the input from the coupler are done at the same time step as the output, then we could use the TFN array and properly increment its values ...
OK, so we could increase the dimension of FLAGS from 1:8 to 1:9 and add a last input type which would be the coupling and process this new input like the other ones. The problem is that this means changing the format of all ww3_shel.inp files ... A NIGHTMARE.
The other options is to do it internally: define the settings for this 8th input type form the information on the 7th output type. So we on,ly read from ww3_shel.inp the FLAGS(1:7) and then we set FLAGS(8) and the TFN(8) from the ODAT ...
But, the forcing from the coupling will also be current and water level, like the J=2 and J=1 input types. So this is not really an extra input, and imagine (let's be crazy) that you only want to couple currents and still read the levels from the usual level.ww3 file.
So it may be better to introduce a bunch of statements with !/COU to manage the different forcing fields separately and set the values of TFN from data received from the coupler.
a coupling solution
Let's try something like this
- we leave the loop DO J=1,8 as it is.
- right before this piece of code
! ! 7.b Run the wave model for the given interval ! TIME0 = TTIME
we insert some coupler calls that resets the value of TTIME to the next coupler time, for example,
!/COU NEXT_TIME = TTIME !/COU IF (MYCOUPLER.EQ.'OASIS') CALL GET_OASIS(NEXT_TIME) !/COU IF (MYCOUPLER.EQ.'ESMF') CALL GET_ESMF(NEXT_TIME) !/COU DTTST = DSEC21 (TTIME , NEXT_TIME ) !/COU IF ( DTTST.GT.0.) TTIME= NEXT_TIME
we may add some further safety checks on NEXT_TIME which should be larger than TIME ... The result should be that the next W3WAVE integration will not carry us further than the next rendez-vous with the coupler. We may also need to check that the changes on TTIME done for J=1 and J=2 (as if the model was reading levels and currents from files) are not making a mess.
This is for the input. For the output, we have to be sure that the parameters to be pushed to the coupler are indeed computed. For some of them (the fluxes like TAUOC, PHIOC ... ) no problem, these are always updated when the model is integrated (although we may prefer to have CUMULATED fluxes and not instataneous values ... , this could require an extra buffer variable that computes the cumulated values inside of W3WAVE, after each call to W3SRCE) . For the others, they will only be computed if the output / coupler step is properly defined. With this logic the "input times" are the master clock for the coupler, and WWATCH will return whatever it will have computed.