Numerical simulation of waves over the ocean.

What is a wave model?

 

·  Waves 

There are many types of waves in the ocean. Probably the first type of wave that comes to mind is a wind wave (wind wave) formed by wind, as shown below. (There are also waves that originate from the depth of the ocean and topography (such as earthquakes), but I will focus on the most common ones here. In terms of wave theory, wind waves are waves that are in the process of developing due to wind. Therefore, the stronger the wind is, the bigger the wind waves become. The reason why experienced surfers and mariners always check the wind at sea is to know the condition of the waves. Waves that have similar characteristics to wind waves at sea are called swells. A swell is a wave when the wind that causes the wind wave stops, or when the wind wave moves into an area where there is no wind. From a wave theory point of view, wind waves are in the process of development, while swells are in the process of decay and progress slowly. One of the characteristics of swell is that the longer the wavelength of the wave, the faster the phase speed. This feature explains why typhoon-derived wind waves generated at low latitudes transform into swells over strong wind areas and reach the coast of Japan.

 

These wind waves and swells together are called waves, and a wave model is a numerical simulation of wind waves and swells based on wind (sea breeze) data. 

·  Wave model

 

When simulating waves, a major problem is the difference in grid spacing between the real space and the model. A grid spacing of several tens of meters is required to represent the shape of a wave where many waves interact with each other. However, when simulating the ocean, it is rare to simulate only a single point. Therefore, it is not practical to analyze each wave individually when simulating a vast ocean. The important point of wave models is that they do not calculate each wave directly, but use wave spectra that apply the properties of waves to infer changes in waves from changes in their statistics. The wave spectrum is expressed as E = (θ, σ), where E is the energy of the wave, θ is the direction of wave propagation (radian), and σ is the frequency of the wave (1/sec). In other words, the wave is represented as a spatial energy distribution by using θ and σ. Therefore, the basic equation of the wave model is the following energy balance equation, which describes the change of energy in the wave spectrum with time. In the following equation, the first term on the left-hand side is the time variation of wave energy, and the second term on the left-hand side, (nabla), indicates the operator that maps the energy value of the vector field to the scalar field we are discussing. The second term on the left-hand side, (nabla), represents the operator that maps the energy value of the vector field to the scalar field under discussion. The term on the right-hand side describes the external forces related to the development, decay, and dissipation of waves, and is commonly referred to as the source function, where Snet means the net amount of external forces that affect the shape of the wave.

  The first is the energy input term Sin due to the ocean wind, the second is the energy dissipation term Sds due to wave dissipation and attenuation such as breaking, and the third is the energy transport term Snl due to nonlinear interaction. The second is the energy dissipation term Sds due to dissipation and attenuation of waves such as breaking. The third term, the energy transport term due to nonlinear interaction, indicates the energy transport related to wave development. The third term, the energy transport term due to nonlinear interaction, describes the energy transport related to wave development and is nonlinear because it describes the interaction between waves. Recent advances in technology have made it possible to handle this nonlinear energy transport term, and more accurate models have been developed.

 

   The wave model that is widely used in recent research is called the third generation model. Specifically, there are WAVEWATCH developed and operated by NOAA (National Oceanic and Atmospheric Administration), SWAN developed and operated by Delft University of Technology in the Netherlands, and WAM developed and operated by BODC (British Oceanographic Data Centre).

The University of Miami Wave Model (UMWM)

  

The wave model presented here is the University of Miami Wave Model (UMWM). This model is a 2.5 generation wave model. As mentioned above, most of the current research is based on third-generation wave models that take into account nonlinear energy interactions, but the UMWM replaces the nonlinear energy interaction term with a simple bulk equation. UMWM replaces the nonlinear energy interaction term with a simple bulk equation, where bulk means "rough" and the constant is an empirical value obtained from past experimental data. This operation reduces the amount of complicated and computationally expensive nonlinear terms and makes the entire model lighter. In fact, the lightness of the model is the most important feature of UMWM, and since the overall source code is short, it is easy to understand what is being done in the model calculation and to improve it. Also, the accuracy of the model results is calculated at a reasonably high level, making the model very easy to handle.

 

·  Installing UMWM

 

We will now install UMWM in the Linux environment and conduct the experiment. First of all, install (compile) according to this manual. The library used for compilation is NetCDF. (For details, see the LOVECLIM article)

After the installation of NetCDF is finished, compile UMWM. Extract the downloaded UMWM gz file and create the UMWM-1.0.0 directory.

 

Move to the source code directory src under UMWM-1.0.0 and start working. For more information about compiling libraries and models, please refer to the LOVECLIM article above.

The first step to compile is to edit the Makefile in /src.

 

# University of Miami Wave Model Makefile

 

# Compiler, flags, and libraty paths:

 

FC         = mpiifort  -save-temps -DMPI

#FCFFLAGS   = -O3 -fastsse

FCFFLAGS   = -g 

CPPFLAGS   =

LD         = $(FC)

LDFLAGS    = $(FCFFLAGS)

NETCDF     = /home/lib/netcdf-4.1.3

NETCDFPATH = -I$(NETCDF)/include -L$(NETCDF)/lib -lnetcdff

NETCDFLIB  = $(NETCDF)/lib/libnetcdff90.a 

 

The important thing here is to specify the compiler and set the NetCDF path.

As described in the UMWM manual, UMWM can perform calculations in an MPI environment, which is an option to shorten calculations by parallelizing the computer (CPU). So, in this case, we use FC = mpiifort to make the compiler ifort. (-DMPI is a definition that says to use MPI when using the compiler, but I needed this definition in my environment...)

Next, specify the path to NetCDF as the directory where you installed it.

 

This completes the editing of the Makefile.

Finally, enter the make command in /src and compile. If there are no errors, you will find the executable file named umwm under /umwm-1.0.0 will be created under /umwm-1.0.0. If you get an error, modify the Makefile and recompile with make clean, make.

 

 Also, there is a directory called /tools, which is used when the output data (grid and bathymetry data) of the domain meteorological model (WRF) is used in UMWM. This time, we used a different method for our experiments, so we will skip the explanation.

 

 

·  Preparing for the UMWM experiment

After successfully completing the compilation, we will start the experiment.

First of all, the essential input data for running UMWM is the sea surface wind data. Other data such as depth data and sea surface current data can be input optionally, but for this experiment, we will use only the sea surface wind data as input data.

We will use data from MSM, the Japan Meteorological Agency's forecasting model, which has a spatial and temporal resolution of 1 hour for every 5 km. MSM is the best model for reproducing past atmospheric data because it consists of forecast values that fill in the gaps between the original data. The data is available here

The file format is NetCDF (.nc), so you can easily create a figure by creating a control file (.ctl). See the LOVECLIM article above for details.

 

First, let's download the MSM u- and v-component wind data (10m) and make a figure. We want the land surface data on the page of the previous link, so we select MSM in the NetCDF data (MSM, RSM) reconstructed around the analysis values, and then select the MSM-S directory on the next page. Under the year directory, the .nc files for each day are grouped together, so download the .nc file for the period of interest. In this explanation, we will focus on TC0914, which traveled over the southeast coast of Japan in 2009, so we set the time period to about five days and downloaded it.

The following animation is a GrADS drawing of the downloaded ocean wind data.

The control files and scripts used at this time are omitted. (The script is explained below.)"

Now that the data is ready, let's start preparing to run UMWM.

Normally, terrain data (grid point spacing) is required to run UMWM, but this time, since we are using MSM data, the grid points are already fixed. Therefore, we need to convert the downloaded data to NetCDF file, which is the input format of UMWM.

 

First of all, as you can see from the section 6.4 Input files in the manual I mentioned earlier, the name of the NetCDF file you want to create must be umwmin_YYYY-MM-DD_hh:mm:ss.nc. Another file called umwm.gridtopo is also required, and as mentioned above, we will add this information to the input file and write it. The following is a program to rewrite the NetCDF file.

 

The format of this program is .gs, which stands for GrADS script. The .gs format can do a lot more than the .gex format executable.

Now let's look at the program. The first thing to note is that the entire program is enclosed in ' ' , which specifies a loop. By enclosing the entire command, recursive processing is performed in the specified range

 

'reinit'

'open dmsm-s.ctl'                 MSMファイルを開く.ctlファイル

'q dims'

tt=14905                          開始時刻の指定

while (tt<15242)                  終了時刻の指定 

'set x 340 441'                   x方向(経度方向)の範囲指定

'set y 140 241'                   y方向(緯度方向)の範囲指定

'set t 'tt                        時間の指定 

say tt

'q time'

mytime = subwrd(result,3)         ここで組み込み関数を使い、q time で表示した時間の情報から mytime, myhour, myday,  

myhour = substr(mytime,1,2)       mymonth, myyear にそれぞれを指定しています。

myday = substr(mytime,4,2)

mymonth = substr(mytime,6,3)

myyear = substr(mytime,9,4)

if (mymonth = 'JAN')              月の表示を文字から数字へと変換します。

 mymonth = '01'

else 

if (mymonth = 'FEB')

 mymonth = '02'

else 

if (mymonth = 'MAR')

 mymonth = '03'

else 

if (mymonth = 'APR')

 mymonth = '04'

else

if (mymonth = 'MAY')

 mymonth = '05'

else 

if (mymonth = 'JUN')

 mymonth = '06'

else 

if (mymonth = 'JUL')

 mymonth = '07'

else 

if (mymonth = 'AUG')

 mymonth = '08'

else 

if (mymonth = 'SEP')

 mymonth = '09'

else 

if (mymonth = 'OCT')

 mymonth = '10'

else 

if (mymonth = 'NOV')

 mymonth = '11'

else (mymonth = 'DEC')

 mymonth = '12'

endif

endif

endif

endif

endif

endif

endif

endif

endif

endif

endif

say myhour

say myday

say mymonth

say myyear

'uc=ABS(u*0)'                                    From here, we will use the sdfwrite command to create a NetCDF file.

'set sdfwrite -flt tmp1.nc'                      uc stands for 'x-component of ocean surface currents', and

'sdfwrite uc'                                    we will fill in the data as 0 since it is not used in this 'sdfwrite uc' experiment.

**************************

'vc=ABS(v*0)'                                    The same goes for vc. It means the velocity component of the sea surface and the unit is [m/s].

'set sdfwrite -flt tmp2.nc'

'sdfwrite vc'

'!ncks -h -A tmp2.nc tmp1.nc'                    ! is a command on a terminal, adding tmp2.nc to tmp1.nc.

**************************

'aa = const(maskout(1,100-psea+sp),0,-u)'         u10 (wind speed at 10m above sea level) is the value we need this time, so we write it out

'u10=u*aa'                                       We are trying to remove the land data here.

'sdfwrite u10 tmp2.nc'                           The sum of psea (sea surface pressure) and sp (surface pressure) is subtracted from 100 to get the negative value.

'!ncks -h -A tmp2.nc tmp1.nc'                    The missing value 0 is given to the land area and 1 to the sea area where it is positive.                   

**************************                       Finally, multiply by the wind speed value u in the MSM to complete.

'aa = const(maskout(1,100-psea+sp),0,-u)'        const(maskout(1,100-psea+sp),0,-u)' Do the same for v10.

'v10=v*aa'

'sdfwrite v10 tmp2.nc'

'!ncks -h -A tmp2.nc tmp1.nc'

**************************

*'aa = const(maskout(1,100-psea+sp),0,-u)'

*'z=100*aa'

*'sdfwrite z tmp2.nc'

*'!ncks -h -A tmp2.nc tmp1.nc'

************************

'!mv tmp1.nc umwmin_'myyear'-'mymonth'-'myday'_'myhour':00:00.nc'  Now we can use the built-in function we used earlier to create the UMWM

*'!rm tmp*.nc'                                                     to the name of the input format.

*say result

tt=tt+1                                                            specify the loop by advancing the time by one.

endwhile 

'disable fwrite'

say 'finished' 

 

If you run this .gs file using the run command on GrADS, the aforementioned UMWM input format file will be output as hourly data.

 

·  UMWM experiment

 

 Now that the input file is ready, let's start the UMWM experiment.

First, move all the input files you have created under /umwm-1.0.0/input.

Next, we will do the initial configuration of UMWM. go back one directory from the input directory and edit main.nml in /namelists under umwm-1.0.0. Open it with a suitable editor.

Basically, it is explained in the comment out, but I will describe the important parts.

 

&DOMAIN

  isGlobal     = .FALSE.               ! Global (.TRUE.) or regional (.FALSE.)

  mm           = 282                   ! Domain size in x           x・y方向(経度・緯度方向)の範囲指定です。    

  nm           = 282                   ! Domain size in y           inputファイルの大きさと合わせます。 

  om           = 37                    ! Number of frequency bins

  pm           = 24                    ! Number of directions

  f1           = 0.0313                ! Lowest frequency bin [Hz]

  f2           = 2.0                   ! Highest frequency bin [Hz]

  startTimeStr = '2009-09-18_00:00:00' ! Simulation start time      実験開始時時刻と終了時刻の指定です。

  stopTimeStr  = '2009-09-23_00:00:00' ! Simulation end time  

  dtg          = 3600                  ! Global (I/O) time step [s]

  restart      = .TRUE.               ! Restart from file          実験を途中から始めるための中間ファイルの作成の有無

/                                                                  の設定。

&PHYSICS

  g           = 9.80665 ! Gravitational acceleration   [m/s^2]

  nu          = 0.9E-6  ! Kinematic viscosity of water [m^2/s]

  sfct        = 0.07    ! Surface tension

  kappa       = 0.4     ! Von Karman constant 

  z           = 10.     ! Height of the input wind speed [m]

  gustiness   = 0.0     ! Random wind gustiness factor (should be between 0 and 0.2)

  dmin        = 10.     ! Depth limiter [m]

  explim      = 0.7     ! Exponent limiter (0.69 ~ 100% growth)

  sin_fac     = 0.11    ! Input factor from following winds

  sin_diss1   = 0.10    ! Damping factor from opposing winds

  sin_diss2   = 0.01    ! Damping factor from swell overrunning wind

  sds_fac     = 42.     ! Breaking dissipation factor

  sds_power   = 2.5     ! Saturation spectrum power

  mss_fac     = 120.    ! Mean-square-slope adjustment to Sds

  snl_fac     = 5.      ! Wave energy downshifting factor

  sdt_fac     = 0.01    ! Dissipation due to turbulence factor

  sbf_fac     = 0.001   ! Bottom friction coefficient [m/s]

  sbp_fac     = 0.001   ! Bottom percolation coefficient [m/s]

/

&GRID

  gridid  = 2        ! Set 1 for Cartesian or 2 for spherical projection

  delx    = 0.0625   ! Grid spacing in x and y in: meters for gridid=1,      MSMデータの格子間隔を記入します。

  dely    = 0.05     ! degrees for gridid=2

  lat_s   = 24.85    ! Latitude of bottom left point; used for gridid=2      緯度・経度の開始点を記入します。

  lon_w   = 132.4375   ! Longitude of bottom left point; used for gridid=2

/

&INPUT_FLAGS        ! Set to 1 for input from file, or 0 for homogenous

                    ! fields input from &INPUT list below;

  dflag   = 0       ! Bathymetry

  excom   = 0       ! Set to 1 to exclude user chosen basins/lakes

  wflag   = 1       ! Winds                                                 inputファイルを使用するので 1。

  pflag   = 0       ! Surface pressure

  tflag   = 0       ! Air temperature

  cflag   = 0       ! Currents

  rflag   = 0       ! Water density

/

&INPUT

  depth   = 3000    ! Water depth [m]                       深度を3000m に固定。

  u0      = 25      ! Wind magnitude [m/s]

  wdir0   = 0.      ! Wind direction [rad]

  sfcp0   = 101325. ! Surface pressure [Pa] 

  temp0   = 300.    ! Surface air temperature [K]

  uc0     = 0.      ! x-component ocean current [m/s]

  vc0     = 0.      ! y-component ocean current [m/s]

  rhow0   = 1025.   ! Water density [kg/m^3]

/

&OUTPUT

  outgrid     = .TRUE.  ! Gridded output 

  outspec     = .FALSE. ! Spectrum output 

  xpl         =  45     ! Grid cell in x for stdout (screen)

  ypl         =   6     ! Grid cell in y for stdout (screen)

 

  restart_out = .TRUE.  ! Restart output

 

Basically, the default state is fine. If you want to change the physical quantity, you can reflect it by changing the item here.

 

Now, in the /umwm-1.0.0 directory, you can type the execution command . /umwm in the /umwm-1.0.0 directory to start the experiment, and now we will talk about the computational complexity of the experiment. This time, as we mentioned in the compilation of UMWM, we will calculate in MPI environment. However, due to the huge amount of computation, it takes a lot of time even if we use the MPI environment (so-called supercomputer), so it will take a lot of time if we don't use the MPI environment.

 

This is why we set the restart option to TRUE as described in main.nml, so that we can start the calculation from the middle.

 

 Since the specifications of supercomputers differ from machine to machine, it is impossible to give a definitive explanation, but basically, you set up the nodes (CPU and memory chunks) and cores (CPU), and then perform the calculation with the submit command.

In my example, I used 2 nodes and 24 cores to run the calculation for 1 hour, because I wanted to shorten the time until the supercomputer started the calculation, and also to fix bugs as soon as they appeared.

 

·  Experimental results of UMWM

If the calculation finishes without any problem, you will find a NetCDF file named umwmout- in the /umwm-1.0.0/output directory. (My personal reminder for this experiment is that the value of output data at the start time of calculation is 0, so we need to deal with it. Now, let's create a .ctl file and draw the figure.

 

The .ctl file is here.

 

Here is the animation we will create.

 

The first thing to explain about the animation is that the color indicates the Significant wave height. Significant wave height is the average of significant waves, and the unit is [m]. Significant wave height is the average of the wave height and period of waves observed during a certain period of time, averaged over 1/3 of the total waves in order from highest to lowest. This is what is often referred to as wave height in weather forecasts.

The animation also shows the sea breeze, which is indicated by the arrows, so you can clearly see that the significant wave height is moving as the typhoon progresses. You can also see that the wind is stronger on the right side of the typhoon, where the wind blows in, and the significant wave height is higher.

 

Lastly, I will explain the script for creating this animation. Here is the script:

  

 As in the case of rewriting the NetCDF file, a loop is used, so it is enclosed in ''. Also, the script is almost the same as when I depicted the MSM data, so please refer to it.

 

'reinit'       

'open umwmout.ctl'       

tt=121                           開始時刻の指定をします。

'set display color white'  

while (tt<181)                   終了時刻の指定をします。

'c'         

'set t 'tt  

'q dims'

say tt

result2 = sublin(result,5)      作成する図に時刻を表示させるため、組み込み変数を使用します。

say result2

mydate = subwrd(result2,6)

say mydate

'set lon 133 150'               経度の領域を指定します。

'set lat 25 39'                 緯度の領域を指定します。

'set grid off'             

'set grads off'            

'set rgb 19 130 0 220'     

'set rgb 20 30 60 255'     

'set rgb 21 0 160 255'     

'set rgb 22 0 200 200'     

'set rgb 23 0 210 140'     

'set rgb 24 0 220 0'       

'set rgb 25 160 230 50'    

'set rgb 26 230 220 50'    

'set rgb 27 230 175 40'    

'set rgb 28 240 130 40'    

'set rgb 29 250 95 50'     

'set rgb 30 250 60 60'     

'set rgb 31 240 0 130'     

'set rgb 32 240 0 200'     

'set rgb 33 240 0 240'     

'set clevs 0 1 2 3 4 5 6 7 8 9 10 11'                    有義波高の色分けの設定です。他にも様々な方法があるので

'set ccols 20 21 22 23 24 25 26 27 28 29 30 31 32'       いろいろと調べてみてください。

'set gxout shaded'                                       色分けの設定です。

'd swh'                                                  有義波高をここで描写しています。

'run cbarn.gs 0.9 1 9.7 4.3'                             カラーバーの表示設定です。

'draw ylab Latitude'

'draw xlab Longitude'

'set ccolor 1'

'set cthick 10'

'set arrscl 1 40' 

'd skip(uwnd,12,12);skip(vwnd,12,12)'                      ここで、矢印の設定をした後、inputした風速を表示しています。

'draw title Swh (color) and Wind speed (arrows) 'mydate''  

'basemap.gs L 0'                                         basemap.gsというスクリプトを使用して陸域をマスクしています。

'printim im-'tt'.png x1000 y800'                         画像の保存をしています。

tt=tt+1                                                  ループの設定をしています。

endwhile    

'!convert -delay 20 im-* swh.gif'                        アニメーションの作成コマンドです。

'!rm im-*'

 

 We use basemap.gs to mask the land area, but the script forthis (color scale, etc.) is available on the net, so you may want to check it out.