module mod_common_wavesice !! contains global wave variables, and the following subroutines !! and functions: !! !! subroutine waves_fsd_reset !! subroutine waveadv_lax(g,ag,dirn,scpy,scpx,dt) !! subroutine waveadv_weno(h,u,v,scuy,scvx,scp2i,scp2,dt) !! subroutine weno3pd_v2(g,sao,u,v,scuy,scvx,scp2i,scp2,dt) !! subroutine dump_netcdf(ncfil,fld_loc,vbl_name) !! subroutine set_inc_waves() !! function theta_dirfrac(th1,th2,mwd) #if defined(WAVES) use mod_xc implicit none !! TODO: change nbdy to margin throughout; !! - need to be careful as things like ficem and hicem are defined with !! nbdy; !! define number of wave frequencies and directions !! TODO define frequencies and directions here !! (currently in mod_waves_init::waves_init) integer, parameter :: n_wavefreq = 25 #if defined (WAVES_ALL_DIRS) integer, parameter :: n_wavdir = 16 !! theta range: !! avgmwd-dth_tot/2.dfloe_miz real,parameter :: dfloe_pack_init = 300.0!unbroken ice initially has this value real,parameter :: dfloe_pack_thresh = 250.0!dfloe can drop due to advection, giving the appearance of breaking private :: theta_dirfrac contains #if defined (WAVES_LAXADV) c Ice advection routine c c subroutine waveadv_lax(g,ag,dirn,scpy,scpx,dt) ! ! --- ------------------------------------------------------------------ ! --- LB 08.06.11 ! --- The Lax scheme for the advection eqn in 2d space ! --- TW 20140717 ! --- Added fix to kill waves on land with LAND_MASK (at the end) ! --- TODO: need to add flux limiting as there is too ! --- much numerical diffusion; ! --- ------------------------------------------------------------------ ! use mod_xc ! implicit none ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & g,ag,scpy,scpx real, dimension(1:idm,1:jdm) :: & gx,gy,gtil real dt,dirn,cd,sd,cd2,sd2 real, parameter :: pi = 3.141592653589793 ! integer i,j,l,im1,im2,ip1,jm1,jm2,jp1 integer itest,jtest itest=327 !376 !290 jtest=132 !49 !30 cd = cos(pi*(dirn-180)/180)!uwave sd = sin(pi*(dirn-180)/180)!vwave cd2 = cd**2 sd2 = sd**2 ! check the value at itest, jtest if ( ((itest.gt.i0) .and. (itest.lt.(i0+idm))) .and. & ((jtest.gt.j0) .and. (jtest.lt.(j0+jdm))) ) then print*, 'LBtest (LAX i): vel, dt, cos, sin, dx, dy', & ag(itest,jtest), dt, ! & cd, sd, scpx(itest,jtest), scpy(itest,jtest) ! print*, 'LBtest (LAX ii): ii, jj, i0, j0', & ii, jj, i0, j0 print*, 'LBtest (LAX i.ii): length(isp;ifp,ilp)', & size(isp), size(ifp,1), size(ilp,1) print*, 'LBtest (LAX i.iii): tests', & scpx(1-nbdy,1-nbdy), ifp(1-nbdy,1), ilp(1-nbdy,1) c print*, 'LBtest (LAX i.iv): size(scpx;scpy)', c & size(scpx,1), size(scpx,2), size(scpy,1), size(scpy,2) end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! LB: isp determines whether the cell is land or not ! ! Calculate derivatives do j=1,jj ! if ( ((itest.gt.i0) .and. (itest.lt.(i0+idm))) .and. ! & (jtest-j0.eq.j) ) then ! print*, 'LBtest (LAX ii.i): ', isp(j), ! & ifp(j,1), ilp(j,1) ! print*, 'LBtest (LAX ii.ii): ', ! & ifp(j,2), ilp(j,2) ! end if c do l=1,isp(j) c do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) do i=1,ii if ((itest-i0.eq.i) .and. (jtest-j0.eq.j)) hen print*, 'LBtest (LAX iii): looping', i, j end if gx(i,j)=sd*ag(i,j)* & (g(i+1,j)-g(i-1,j))/(scpx(i,j)+ & 0.5*(scpx(i+1,j)+scpx(i-1,j))) end do c end do end do c if ( ((itest.gt.i0) .and. (itest.lt.(i0+idm))) .and. c & ((jtest.gt.j0) .and. (jtest.lt.(j0+jdm))) ) then c print*, 'LBtest: -> gx = ', gx(itest,jtest) c end if do j=1,jj c do l=1,isp(j+nbdy) c do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) do i=1,ii gy(i,j)=cd*ag(i,j)* & (g(i,j+1)-g(i,j-1))/(scpy(i,j)+ & 0.5*(scpy(i,j+1)+scpy(i,j-1))) end do c end do end do c if ( ((itest.gt.i0) .and. (itest.lt.(i0+idm))) .and. c & ((jtest.gt.j0) .and. (jtest.lt.(j0+jdm))) ) then c print*, 'LBtest: -> gy = ', gy(itest,jtest) c end if ! Spatial average do j=1,jj c do l=1,isp(j) c do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) do i=1,ii if ((itest-i0.eq.i) .and. (jtest-j0.eq.j)) then print*, 'LBtest (LAX iv): gi', & sd2, g(i+1,j), g(i-1,j) print*, 'LBtest (LAX v): gj', & cd2, g(i,j+1), g(i,j-1) end if gtil(i,j)= 0.5*( & sd2*(g(i+1,j)+g(i-1,j)) + & cd2*(g(i,j+1)+g(i,j-1)) ) end do c end do end do c if ( ((itest.gt.i0) .and. (itest.lt.(i0+idm))) .and. c & ((jtest.gt.j0) .and. (jtest.lt.(j0+jdm))) ) then c print*, 'LBtest: -> gtil = ', gtil(itest,jtest) c end if ! Advect do j=1,jj c do l=1,isp(j) c do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) do i=1,ii g(i,j)=gtil(i,j) - dt*(gx(i,j)+gy(i,j)) end do c end do end do where (LAND_MASK==1.0) g = 0.0 return end subroutine waveadv_lax #else !!old flag (WAVES_ADV_WENO) - now this is the default c Ice advection routine - adapted for waves c subroutine waveadv_weno(h,u,v,scuy,scvx,scp2i,scp2,dt) ! ! --- ------------------------------------------------------------------ ! --- Advection is done with flux limited 3rd order WENO in space and ! --- 2nd order Runge-Kutta in time ! --- ------------------------------------------------------------------ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! modified [TW 28.1.2014] for wave advection !! - iceadv builds up ice (& thus waves) at coast; !! - waveadv_weno lets waves go into/out of land, !! then cancels the waves on land at end of routine !! (with LAND_MASK); !! - just adds the wave cancellation to Luke's original fix !! for waves on boundaries; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! ARGUMENTS: !! in/out: h is qty to be advected; !! !! in: u,v,dt are speed components and time step; !! in: scuy,scvx is mesh size at u (v) points in y (x) direction !! - see common_blocks.h; !! in: scp2, scp2i are grid box area at p points, and its inverse; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use mod_xc ! implicit none ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & h,u,v,scuy,scvx,scp2i,scp2 real dt ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: sao,hp real dtm integer i,j,l sao(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)=0 ! ! --- Use a modified time step since velocities are in m/s while scale ! --- factors are in cm ! dtm=dt*1.e2 dtm=dt ! --- Prediction step call weno3pd_v2(h,sao,u,v,scuy,scvx,scp2i,scp2,dtm) !$OMP PARALLEL DO margin=nbdy do j=1-margin,jj+margin ! do l=1,isp(j) ! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do i=1-margin,ii+margin!![TW 28.1.2014] hp(i,j)=h(i,j)+dtm*sao(i,j) enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO !TW 20151021: do tiling here !TODO this is the easy but slow way - faster thing would be to !change the loops in weno3pd_v2 !TODO is this necessary - if nbdy=6, maybe errors from boundary ! can't creep in in 1 time step? call xctilr(hp (1-nbdy,1-nbdy),1,1, 6,6, halo_ps) ! ! --- Correction step call weno3pd_v2(hp,sao,u,v,scuy,scvx,scp2i,scp2,dtm) !$OMP PARALLEL DO do j=1-margin,jj+margin ! do l=1,isp(j) ! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do i=1-margin,ii+margin!![TW 28.1.2014] if (LAND_MASK(i,j).eq.0.0) then !!if not on land keep as it is; h(i,j) = .5*(h(i,j)+hp(i,j)+dtm*sao(i,j)) else !!if on land set quantitity (waves) to zero h(i,j) = 0.0 end if enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO return end subroutine waveadv_weno subroutine weno3pd_v2(g,sao,u,v,scuy,scvx,scp2i,scp2,dt) ! ! --- ------------------------------------------------------------------ ! --- By a weighted essentially non-oscillatory scheme with up to 3th ! --- order accuracy, obtain the spatial advective operator of a ! --- 2-dimensional field defined at the scalar points of a C-grid. The ! --- fluxes are limited to make the scheme positive definite. ! --- Advective velocities in the i- and j-direction are defined at u- ! --- and v-points, respectively. ! --- ------------------------------------------------------------------ !! TW 20140128: changed to allow waves onto land, !! and kill them at the end of the routine with the LAND_MASK !! - iceadv sets u,v=0 on land so ice accumulates at coasts ! use mod_xc ! implicit none ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & g,sao,u,v,scuy,scvx,scp2i,scp2 real dt ! real cq00,cq01,cq10,cq11,ca0,ca1,eps parameter (cq00=-1./2.,cq01= 3./2., & & cq10= 1./2.,cq11= 1./2., & & ca0=1./3.,ca1=2./3., & & eps=1.e-12) ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & ful,fuh,fvl,fvh,gt real q0,q1,a0,a1,q integer i,j,l,im1,im2,ip1,jm1,jm2,jp1 ! ! --- Compute grid cell boundary fluxes. Split in a low order flux ! --- (donor cell) and a high order correction flux. ! !$OMP PARALLEL DO do j=0,jj+2 do i=0,ii+2 ful(i,j)=0. fuh(i,j)=0. fvl(i,j)=0. fvh(i,j)=0. enddo enddo !$OMP END PARALLEL DO ! ! call xctilr(g, 1,1, 3,3, halo_ps) ! !$OMP PARALLEL DO PRIVATE(im1,im2,q0,q1,a0,a1,ip1) do j=0,jj+1 ! do l=1,isu(j) ! do i=max(0,ifu(j,l)),min(ii+2,ilu(j,l)) do i=0,ii+2!![TW 28.1.2014] im1=i-1 ! if (u(i,j).gt.0.) then !iu is a water mask (water at point and to left of point); !for waves we make it 1 everywhere and mask waves that go on land later; !im2 = im1-iu(im1,j) im2=im1-1 ! q0=cq00*g(im2,j)+cq01*g(im1,j) q1=cq10*g(im1,j)+cq11*g(i ,j) ! a0=ca0 a1=ca1*(abs(g(im2,j)-g(im1,j))+eps) & & /(abs(g(im1,j)-g(i ,j))+eps) ! ful(i,j)=u(i,j)*g(im1,j)*scuy(i,j) ! else !ip1=i+iu(i+1,j) ip1=i+1 ! q0=cq11*g(im1,j)+cq10*g(i ,j) q1=cq01*g(i ,j)+cq00*g(ip1,j) ! a0=ca1 a1=ca0*(abs(g(im1,j)-g(i ,j))+eps) & & /(abs(g(i ,j)-g(ip1,j))+eps) ! ful(i,j)=u(i,j)*g(i ,j)*scuy(i,j) ! endif ! fuh(i,j)=u(i,j)*(a0*q0+a1*q1)/(a0+a1)*scuy(i,j)-ful(i,j) ! enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO ! !$OMP PARALLEL DO PRIVATE(jm1,q0,q1,a0,a1,jm2,jp1) do j=0,jj+2 jm1=j-1 ! do l=1,isv(j) ! do i=max(0,ifv(j,l)),min(ii+1,ilv(j,l)) do i=0,ii+1!![TW 28.1.2014] ! if (v(i,j).gt.0.) then !iv is a water mask (water at point and beneath point); !for waves we make it 1 everywhere and mask waves that go on land later; !jm2=jm1-iv(i,jm1) jm2=jm1-1 ! q0=cq00*g(i,jm2)+cq01*g(i,jm1) q1=cq10*g(i,jm1)+cq11*g(i,j ) ! a0=ca0 a1=ca1*(abs(g(i,jm2)-g(i,jm1))+eps) & & /(abs(g(i,jm1)-g(i,j ))+eps) ! fvl(i,j)=v(i,j)*g(i,jm1)*scvx(i,j) ! else !jp1=j+iv(i,j+1) jp1=j+1 ! q0=cq11*g(i,jm1)+cq10*g(i,j ) q1=cq01*g(i,j )+cq00*g(i,jp1) ! a0=ca1 a1=ca0*(abs(g(i,jm1)-g(i,j ))+eps) & & /(abs(g(i,j )-g(i,jp1))+eps) ! fvl(i,j)=v(i,j)*g(i,j )*scvx(i,j) ! endif ! fvh(i,j)=v(i,j)*(a0*q0+a1*q1)/(a0+a1)*scvx(i,j)-fvl(i,j) ! enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO ! ! --- Update field with low order fluxes. !$OMP PARALLEL DO do j=0,jj+1 ! do l=1,isp(j) ! do i=max(0,ifp(j,l)),min(ii+1,ilp(j,l)) do i=0,ii+1!![TW 28.1.2014] gt(i,j)=g(i,j)-dt*(ful(i+1,j)-ful(i,j) & & +fvl(i,j+1)-fvl(i,j))*scp2i(i,j) enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO ! ! --- Obtain fluxes with limited high order correction fluxes. q=.25/dt !$OMP PARALLEL DO do j=1,jj ! do l=1,isu(j) ! do i=max(1,ifu(j,l)),min(ii+1,ilu(j,l)) do i=1,ii+1!![TW 28.1.2014] fuh(i,j)=ful(i,j)+max(-q*gt(i ,j)*scp2(i ,j), & & min( q*gt(i-1,j)*scp2(i-1,j),fuh(i,j))) enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO !$OMP PARALLEL DO do j=1,jj+1 ! do l=1,isv(j) ! do i=max(1,ifv(j,l)),min(ii,ilv(j,l)) do i=1,ii!![TW 28.1.2014] fvh(i,j)=fvl(i,j)+max(-q*gt(i,j )*scp2(i,j ), & & min( q*gt(i,j-1)*scp2(i,j-1),fvh(i,j))) enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO ! ! --- Compute the spatial advective operator. !$OMP PARALLEL DO do j=1,jj ! do l=1,isp(j) ! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) do i=1,ii!![TW 28.1.2014] sao(i,j)=-(fuh(i+1,j)-fuh(i,j)+fvh(i,j+1)-fvh(i,j))*scp2i(i,j) enddo ! enddo ! enddo enddo !$OMP END PARALLEL DO ! return end subroutine weno3pd_v2 #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Routine to make it easier to dump netcdf files in an easily !! visualisable way subroutine dump_netcdf(ncfil,fld_loc,vbl_name) use mod_xc use netcdf use mod_hycom_nersc, only: ncdraft,nc_put_att implicit none real, intent(in), & & dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: fld_loc real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: tmp_ij character(len=*), intent(in) :: ncfil !name of netcdf file character(len=*), intent(in) :: vbl_name !name of variable real, dimension(itdm,jtdm) :: fld real, parameter :: netcdf_landval = -32767. real :: tmpval real,save :: tmp_min,tmp_max !! integer :: i,j logical :: crit tmp_ij = fld_loc tmp_min = minval(tmp_ij) call xcminr_0o(tmp_min) tmp_max = maxval(tmp_ij) call xcmaxr_0o(tmp_max) if (mnproc==1) then print*,'Dumping variable to netcdf file:'//trim(vbl_name) print*,trim(ncfil) print*,'Variable min/max:' print*,tmp_min,tmp_max print*,' ' end if do j=1,jj do i=1,ii crit = (.not.(tmp_ij(i,j).ge.0.)) & .and.(.not.(tmp_ij(i,j).le.0.)) if (crit) then print*,trim(vbl_name)//' is NaN at',i+i0,j+j0 end if end do end do where (ip.eq.0) tmp_ij = netcdf_landval call xcaget(fld,tmp_ij,0) call ncdraft(trim(ncfil),fld,trim(vbl_name),'') !!add "missing_value" attribute call nc_put_att(trim(ncfil),trim(vbl_name) & & ,'missing_value',netcdf_landval) call nc_put_att(trim(ncfil),trim(vbl_name) & & ,'_FillValue',netcdf_landval) end subroutine dump_netcdf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine set_inc_waves() implicit none include 'common_blocks.h' !itest,jtest integer, parameter :: nw = n_wavefreq integer, parameter :: ndir = n_wavdir !! real,parameter :: gravity = 9.81 real :: period & ,chi,theta_fac,exB,fac & ,aa,bb,dtheta real :: mom0_test integer :: i,j,w,wth logical :: set_wave_field #if defined (WAVES_TEST_UNIFWAVES) !new flag [TW 13.5.2011] ! uniform waves across water regions ! (instead of getting waves from WAM); real, parameter :: swh_unif = 3.0 ! sig wave height real, parameter :: mwp_unif = 8.5 !mean wave period real :: mwd_unif !mean wave direction logical, parameter :: shrink_lonlatrng = .true. !only set spectrum in a specified range of lon & lat logical :: inlonrng,inlatrng real, parameter :: lon_min = -5.0 real, parameter :: lon_max = 15.0 real, parameter :: lat_min = 73.0 real, parameter :: lat_max = 76.0 #endif #if defined (WAVES_TIME_INTERP) logical :: set_wave_field2 #endif aa = 8.1e-3*gravity!!for Bretschneider spectrum bb = 1.25!!for Bretschneider spectrum if (mnproc==1)print*, 'TW: setting incident wave field' WAVE_MASK = 0.0 #if defined (WAVES_TIME_INTERP) WAVE_MASK2 = 0.0 #endif do j=1,jj do i=1,ii set_wave_field = (mwp(i,j) .gt. 0.0)!some waves defined there by input wave model & .and.(ICE_MASK(i,j).eq.0.0)!no ice there & .and.(LAND_MASK(i,j).eq.0.0)!no land there #if defined (WAVES_TEST_UNIFWAVES) if (set_wave_field.and.shrink_lonlatrng) then !set wave field only within a specified lat & lon range; inlonrng = ( plon(i,j).gt.lon_min ) inlonrng = ( inlonrng.and.(plon(i,j).lt.lon_max) ) !! inlatrng = ( plat(i,j).gt.lat_min ) inlatrng = ( inlatrng.and.(plat(i,j).lt.lat_max) ) !! set_wave_field = (inlonrng.and.inlatrng) end if #endif if ((i.eq.itest).and.(j.eq.jtest)) then print*,'TW wave output check 1' print*,'overwrite?',set_wave_field print*,'Hs ',swh_out(i,j),swh(i,j) print*,'Tp ',mwp_out(i,j),mwp(i,j) print*,'mwd',mwd_out(i,j),mwd(i,j) end if if (set_wave_field) then #if defined (WAVES_TEST_UNIFWAVES) swh(i,j) = swh_unif mwp(i,j) = mwp_unif mwd(i,j) = mwd_unif #endif WAVE_MASK(i,j) = 1.0 swh_out (i,j) = swh(i,j) mwp_out (i,j) = mwp(i,j) mwd_out (i,j) = mwd(i,j) end if if ((i.eq.itest).and.(j.eq.jtest)) then print*,'TW wave output check' print*,'i,j ',i+i0,j+j0 print*,'Hs ',swh_out(i,j) print*,'Tp ',mwp_out(i,j) print*,'mwd',mwd_out(i,j) print*,'WAVE_MASK',WAVE_MASK(i,j) end if #if defined (WAVES_TIME_INTERP) set_wave_field2 = (mwp2(i,j) .gt. 0.0)!waves present in future reading & .and.(set_wave_field) !waves also in the current reading if (set_wave_field2) WAVE_MASK2(i,j) = 1.0 #endif end do!i end do!j dtheta = abs(wavdir(2)-wavdir(1)) mom0_test = 0. do w=1,nw do wth=1,ndir period = wave_period(w) do j=1,jj do i=1,ii if (WAVE_MASK(i,j).eq.1.0) then !theta_fac = 2/pi*cos(chi)**2!!NB do not integrate yet!! ! integrate, then divide by dtheta to give a more ! accurate representation of energy in directional bin theta_fac = theta_dirfrac(wavdir(wth)-dtheta/2., & wavdir(wth)+dtheta/2.,mwd(i,j))/wt_theta(wth) fac = bb*period**5/(8*pi*mwp(i,j)**4) exB = exp(-bb*(period/mwp(i,j))**4) !! sdf_out(i,j,wth,w) = fac*swh(i,j)**2* & exB*theta_fac if ((i.eq.itest).and.(j.eq.jtest)) then mom0_test = mom0_test & +wt_omega(w)*wt_theta(wth)*sdf_out(i,j,wth,w) end if if (mwp(i,j).le.0.0) then print*,'mwp<=0 inside wavemask' print*,'i,j',i+i0,j+j0,mwp(i,j),WAVE_MASK(i,j) & & ,LAND_MASK(i,j) !call xcstop('(bad mwp in set_inc_waves)') end if end if #if defined (WAVES_TIME_INTERP) if (WAVE_MASK2(i,j).eq.1.0) then ! integrate, then divide by dtheta to give a more ! accurate representation of energy in directional bin theta_fac = theta_dirfrac(wavdir(wth)-dtheta/2., & wavdir(wth)+dtheta/2.,mwd(i,j))/wt_theta(wth) fac = bb*period**5/(8*pi*mwp2(i,j)**4) exB = exp(-bb*(period/mwp2(i,j))**4) !! sdf_dir2(i,j,wth,w) = fac*swh2(i,j)**2* & exB*theta_fac end if #endif end do!i end do!j end do!wth end do!w if ((itest.gt.0).and.(jtest.gt.0)) then i = itest j = jtest w = nint(nw/2.) print*,'TW sdf_out (set_inc_waves)' print*,'i,j,w: ',i+i0,j+j0,w do wth=1,ndir print*,'wth, th, S: ',wth,wavdir(wth),sdf_out(i,j,wth,w) end do!wth mom0_test = 4*sqrt(mom0_test)!convert integral to Hs print*,'TW: test integral of inc S: swh=',mom0_test,swh(i,j) end if end subroutine set_inc_waves !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function theta_dirfrac(th1,th2,mwd) !! chi=pi/180*(theta-mwd) !! (convert to radians and centre about mwd, the mean wave direction) !! chi1=pi/180*(th1-mwd) !! chi2=pi/180*(th2-mwd) !! D(chi)=2/pi*cos^2(chi) if chi\in[-pi/2,pi/2] !! D(chi)=0 if chi\in[-pi/2,pi/2] (no backwards waves) !! theta_dirfrac = \int_{chi1}^{chi2}D(chi)dchi !! = 1 if -chi1=chi2=pi/2 !! = the fraction of the energy going in the directions in [chi1,chi2] implicit none real, intent(in) :: th1,th2,mwd real :: chi1,chi2 real, parameter :: pi = 3.1415927 chi1 = pi*(th1-mwd)/180 !chi1 = .5*(-pi/2+chi1)+.5*abs(-pi/2-chi1)!max(-pi/2,chi1) chi1 = min(pi/2,max(-pi/2,chi1))!-pi/2<=chi1<=pi/2 chi2 = pi*(th2-mwd)/180 !chi2 = .5*(pi/2+chi2)-.5*abs(pi/2-chi2)!min(pi/2,chi2) chi2 = min(pi/2,max(-pi/2,chi2))!-pi/2<=chi2<=pi/2 theta_dirfrac = 2*(chi2-chi1)+sin(2*chi2)-sin(2*chi1) theta_dirfrac = theta_dirfrac/2/pi end function theta_dirfrac #endif/* WAVES */ end module mod_common_wavesice ! ----------------------------------------------------------------------------------- ! REVISION HISTORY ! ----------------------------------------------------------------------------------- ! 20140717: Obsolete options removed at svn Revision 328 ! WAVES_FREQSPEC - only frequency spectrum considered (>1 direction not possible) ! WAVES_ENSEMBLE - approx multiple directions by taking mean ! of single direction runs ! ! Also removed WGM method (Dumont et al, 2011) ! - old default option ! & WAVES_NOSTRESSCON (WGM with only strain breakage) ! ! WAVES_DIRSPEC option is now default (multiple directions, including 1) ! ! 20140717: Advection options changed at svn Revision 329 ! WAVES_LAXADV corrected for land issues, but still needs flux limiting ! WAVES_ADV_WENO waveadv_weno is improved version of iceadv ! - corrects for land in a better (and simpler) way ! - pretends no land initially, then kills waves on land later ! (with LAND_MASK) ! ! 20151108: ! New subroutines at svn Revision 400: ! set_inc_waves, dump_netcdf ! Call xctilr between prediction and correction step ! - is this necessary - if nbdy=6, maybe errors from boundary ! can't creep in in 1 time step? ! ! -----------------------------------------------------------------------------------