module mod_wavesice #if defined(WAVES) use mod_xc use mod_common_wavesice, only: n_wavefreq,n_wavdir & ,taux_wav,tauy_wav & ,taux_wav_ice,tauy_wav_ice & ,taux_wav_ocn,tauy_wav_ocn & ,wt_omega,wt_theta implicit none private real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy & & ,n_wavefreq) :: wlng_ice real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy & ,n_wavdir) :: uwave,vwave !!TODO move uwave,vwave setup to waves_init? !!TODO cover 360 deg of direction on model grid & convert to !! lon-lat afterwards? real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & dmax,dave_store,Hs,Tp,Mdir,MASKS real, dimension(n_wavefreq) :: omega,wlng,wspd character(len=12) :: dtime_waves_human!!yyyy_ddd.xxx, where dd starts at 0 integer :: w_test,wth_test public :: waves_advect contains subroutine waves_advect(dtime) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! **Routine for estimating wave field and floe size distribution ! in the ice. It advects & attenuates an incident wave ! field, and estimates floe breaking as the waves travel in; ! **Some updates from Dumont et al (2011) to Williams et al (under ! revision, 2012) !! ! Outputs: ! **'dfloe' - a global variable, which has the maximum floe ! size in it ! - this info is stored as local variable 'dmax' during the routine, ! & it is updated after the routine has finished; ! NB dfloe is defined in mod_common_ice.F; ! **swh_out - significant wave height); ! **mwp_out - period between wave crests (or downward zeros in ! wave elevation); ! - these are dumped into 'wavesice_outfld[dayno].nc', with some ! other quantities; ! TODO: these wave outputs are not quite right when not all ! directions are considered !! ! Attenuation model is changed by the 'ACoption' variable; ! - can use Kohout & Meylan (2008) (ACoption=0), ! - or Bennetts & Squire (2012) (ACoption=1 or 2) ! ('2' has extra Robinson-Palmer damping), ! - or of Williams et al (under revision, 2012) (ACoption=4) ! - latter is an example of a bad attenuation model; ! - ACoption=3 gives the no-draft version of '1'; ! - ACoption=5 gives a different interpolation for the '0' model; ! - models are in the module 'mod_attenuation_ice.F' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !FLAGS FOR THIS ROUTINE: !NB NEED 'WAVES' flag for any of them to work; !Best at the moment: ! -DWAVES ! Flag for the best wave model choice ! (eg -DWAVES_WW3 or -DWAVES_WAMNSEA) ! for the required time period or region ! (default is low-res global WAM); ! ! Model uses full directional spectrum, with number of directions set by ! 'ndir' variable (Default - old WAVES_DIRSPEC flag); ! Single direction options deleted ! (WAVES_FREQSPEC, WAVES_ENSEMBLE, WAVES_NOSTRESSCON) !! !Secondary options: !1. WAVES_TIME_INTERP (old WAVES_SAVE_ICEWAVES2 flag - before 20140718): ! Fastest option? - inside the wave mask interpolate in time ! between wave readings ! - thus there is no need to advect tiles that are completely ! inside the wave mask; ! Waves outside the ice mask are saved till next call to waves_advect; ! Ice conditions (dfloe=max floe size) are exported to mod_common_wavesice.F ! and advected; !2. Default (old WAVES_SAVE_ICEWAVES flag - before 20140718): ! Similar to (1), but no time interpolation between readings; ! Waves outside the wave mask are saved till next call to waves_advect; ! Ice conditions (dfloe=max floe size) are exported to mod_common_wavesice.F ! and advected; ! !Extra options: !WAVES_TESTADV flag -> check advection is working by dumping !WAVES_TEST_READWAVES flag -> check waves are read in correctly ! (into eg 'waves_ww3[dayno].nc'); !WAVES_TEST_UNIFWAVES flag -> uniform wave field across open water ! - can also restrict lon & lat range of waves by setting ! 'shrink_lonlatrng' to '.true.' !WAVES_SKIP flag -> can speed up model if necessary ! by taking every 2nd, 4th or 8th wave reading; ! !WAVES_LAXADV flag -> uses a Lax advection scheme instead of WENO; !-TODO: needs flux limiting to prevent numerical damping ! !WAVES_SAMESPEED flag -> if this isn't here the waves all travel at !the same speed; !-TODO: thin plate ice dispersion relation - use correct wavelength/group velocity; !-TODO: Mindlin plate ice dispersion relation? ! - stops the group velocity becoming infinite for small periods; ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c --- TODO: put in Vernon's dependence on brine vol fraction use mod_xc use mod_hycom_nersc use mod_common_ice use netcdf use mod_attenuation_ice!module-> atten_coeff,damping,ice use mod_readwaves, only: waves_timeres use mod_common_wavesice!swh, mwp etc in there now use mod_wim_prams implicit none include "common_blocks.h" !!******************************************************************** !!** DECLARE VARIABLES: ** !!******************************************************************** real*8 , intent(in) :: dtime !logical :: init_wave !!no of frequencies, from n_wavefreq in mod_common_wavesice.F integer :: nw integer :: ndir !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !nwavesteps not related to 6h interval at the moment !(ie nwavesteps*twave.ne.3600s) !(Dany used 2000 & 50s); integer :: nwavesteps ! Number of wave steps (def'd below) !#if ! defined (WAVES_LAX) ! real, parameter :: twave = 50.0 ! Wave time step (s) - max (CFL)~70s !#else real, parameter :: CFL = .70!0.95 !Courant number (cf. stability) ! real, parameter :: CFL = .27 real :: twave ! real :: cd!= cos(pi*(avg_mwd-180)/180) ! real :: sd!= sin(pi*(avg_mwd-180)/180) ! real :: xmax,ymax real :: period,dx_cfl,amax,dxmin real :: tmin,tmax,alp_nd,alp_dim,dmpg,tmp1,tmp2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & tmp2d,uw,vw !!NB to get wave group velocities need to rotate directions from !'l2m' - lon/lat to model real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,n_wavefreq) :: & ag_eff,ap_eff logical :: set_wave_field real :: alpha,dw real :: dave integer :: i,j,m,n,w integer, dimension(2) :: dumloc ! for testing wave advection; integer :: iyear, ihour, iday character(len=4) :: cyy character(len=3) :: cdd,cdfrac character(len=10) :: vbl_name real, dimension(itdm,jtdm) :: fld #if defined (WAVES_TESTADV) !! if flag defined check advection at regular time-steps; !character(len=10) :: ctt,wstr,varstr integer :: nwrite_testadv, jump_testadv logical :: do_testadv #endif !diagnostic text file character(len=80) :: txtfil integer,parameter :: txtnop = 13 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy & ,n_wavdir,n_wavefreq) :: sdf_dir real :: adv_dir,dtheta,chi real :: dir,ddir integer :: wth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy, & & n_wavefreq) :: alp_nd_all, dmpg_all real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & m0_disp_2d,m2_disp_2d,m0_strain_2d,mth_disp_2d & ,Sfreq,atten_dim,taux_om,tauy_om real :: Tcrest, wlng_crest, & Pcrit,Pstrain real :: mth_test,m0_test,m2_test,veps_test integer :: jcrest real :: om,om1,lam1,lam2,ccrest ! expected no of waves in timestep twave, ! & expected period and wavelength; ! Pcrit=e^{-1} ! Pstrain has to exceed this for breaking to happen ! (prob that crit strain will be exceeded) ! NB not including cavitation at the moment ! as unsure how to implement in a spectral setting ! (non-linear) and also not sure that long waves would produce ! as much stress as the formula predicts; ! NB are still including stress criterion via Hooke's law; !real :: epsc, epsc2 !stress limit & Hooke's law ->epsc2, !epsc=min(epsc1,epsc2) is ultimate strain limit; real :: kice !wave-number in ice; !other spectral variables; real :: theta_fac,avg_mwd real, dimension(n_wavefreq) :: sdf_freq_ij,th_sdf_freq_ij real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: sdf_ij !radial freq, quadrature weights for doing !\int_0^\infty[S(\omega)]d\omega (Simpson's rule) real :: ommin,ommax!limits of radial freq logical :: BC!.true.->ice breaks logical :: IS_LAND!.true.->if waves are going into land ! need to set waves to zero #if defined (WAVES_TIME_INTERP) real,parameter :: t_interp_period = 30*60!reset inc waves every 30mins real :: c_interp_waves,t_interp_waves logical :: DO_INTERP_WAVES #endif !!choose values to go in MASKS (mainly to make the ncview colours in the mask nice :) ) real, parameter :: wave_val = 1.0!wave mask defined and needs to be advected real, parameter :: wave_val2 = 2.0!wave mask at 1st and last end time (-DWAVES_TIME_INTERP) real, parameter :: ice_val = 3.0!ice mask defined real, parameter :: miz_val = 4.0!ice mask defined real, parameter :: land_val = 5.0!land mask defined !! real :: radinv !pi,radian defined in blkdat.F for common_blocks.h !radian = pi/180. : degrees -> radians radinv = 180./pi !radians -> degrees !!******************************************************************** !!** INITIALISE VARIABLES: ** !!******************************************************************** ! if (mnproc==1) then ! print*, 'TW: wa start time: ',dtime ! print*, 'TW: dxmin: ',dxmin,minval(scpx) ! end if nw = n_wavefreq ndir = n_wavdir if (mnproc.eq.1) then print*, 'TW: in waves_advect subroutine:' avg_mwd = wavdir(1+nint(n_wavdir/2.0)) print*, 'avg_mwd, nw, ndir:',avg_mwd & & ,nw,ndir end if !!time of call -> names of output files call forday(dtime,yrflag, iyear,iday,ihour) write(cyy ,'(i4.4)') iyear write(cdd,'(i3.3)') iday-1!!start at iday=0 write(cdfrac,'(i3.3)') int(1.0e3*(ihour/24.0)) dtime_waves_human = cyy//'_'//cdd//'.'//cdfrac !! INITIALIZATION if (mnproc==1) then print*,'' print*,'******************************************************' print*,'TW: parameters for! WIM !?!' print*,'ACoption: ! ',ACoption print*,'AGoption: ! ',AGoption print*,'Brine volume fraction: ',brine_vol print*,'Flexural strength (Pa): ',sigc print*,'Breaking strain: ',epsc print*,'Young^s modulus (Pa): ',young print*,'Damping (Pa.s/m): ',visc_RP print*,'Number of wave directions: ',ndir print*,'******************************************************' print*,'' ! !extra tests of functions: ! print*,'TW fss: f,xi,dmin,dmax',fragility,xi,dfloe_min,150. ! call floe_scaling_smooth(dave,fragility,xi ! & ,dfloe_min,150.,1) ! print*,'TW fss test: mom=1',dave ! call floe_scaling_smooth(dave,fragility,xi ! & ,dfloe_min,150.,2) ! print*,'TW fss test: mom=2',dave end if alp_nd = 0.0 alp_dim = 0.0 dmpg = 0.0 tmp1 = 0.0 tmp2 = 0.0 taux_wav = 0. tauy_wav = 0. wspd = 0.0 wlng = 0.0 ag_eff = 0.0 ap_eff = 0.0 !! dave = dfloe_pack_init dave_store = dfloe_pack_init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MASKS = 0.0 iceh = ICE_MASK*iceh Hs = swh_out Tp = mwp_out Mdir = mwd_out sdf_dir = sdf_out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!check dfloe is within limits !!(it can change with advection) !!TODO check if this is affecting MIZ rheology (thickness pile-up !! at MIZ boundary) where ((LAND_MASK(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy) & +ICE_MASK(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy)).eq.0.0) & dfloe = 0.0!!zero in the ocean/on land where ((dfloe(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).gt.dfloe_pack_thresh) & .and.(LAND_MASK(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).eq.0.0)) & dfloe = dfloe_pack_init where ((dfloe(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).lt.dfloe_min) & .and.(ICE_MASK(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).eq.1.0)) & dfloe = dfloe_min!can drop below dfloe_min due to advection dmax = dfloe !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! omega = 2*pi/wave_period ommin = minval(omega) ommax = maxval(omega) tmin = minval(wave_period) tmax = maxval(wave_period) dw = (ommax-ommin)/(nw-1.0) if (mnproc==1) then print*, 'TWtest: tmax,tmin:' print*, ' ',tmax,tmin end if BC = .false.!breaking criterion value IS_LAND = .false.!check if land is present alp_nd_all = 0.0 dmpg_all = 0.0 wlng_ice = 0.0 !m0_disp = 0.0 !m2_disp = 0.0 !m0_strain = 0.0 !mth_disp = 0.0 Tcrest = 0.0 wlng_crest = 0.0 Pcrit = exp(-1.0)!!this is so E_s>sqrt(2)*\eps_c, to agree !!with monochromatic wave case; Pstrain = 0.0 kice = 0.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! theta_fac = 0.0 dtheta = abs(wavdir(2)-wavdir(1)) adv_dir = 0.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (mnproc==1) then print*, 'TWtest: dw,ommax,ommin,nw:' print*, dw,ommax,ommin,nw end if do w = 1,nw period = 2*pi/omega(w) wlng(w) = gravity*period**2/(2*pi)!!water wavelength !!(ice defined below if AGoption==1 or 2) #if defined (WAVES_SAMESPEED) ! Have waves travel at same speed; wspd(w) = 0.5*sqrt(gravity*wlng(1)/(2*pi)) ! Group speed of fastest wave #else ! Have wave dispersion; wspd(w) = 0.5*sqrt(gravity*wlng(w)/(2*pi)) ! Group speed #endif !!TODO: maybe have ice wave speed/wavelength? ag_eff(:,:,w) = wspd(w) !!group vel ap_eff(:,:,w) = 2*wspd(w)!!phase vel end do ! end spectral loop (w) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!Rotate group velocity vectors from 'l2m' - lon/lat to model; !!Set keeppoint to .true. use only i,j points; if (mnproc==1)print*, & 'TW: rotating directions from lon-lat to model' do wth=1,ndir !! wavdir is waves-from !! convention is 0 from north, -90 from west !! this should -> -pi/2 , 0 adv_dir = -radian*(90.0+wavdir(wth)) uwave(:,:,wth) = cos(adv_dir) vwave(:,:,wth) = sin(adv_dir) call rotate2(uwave(1-nbdy,1-nbdy,wth),vwave(1-nbdy,1-nbdy,wth), & plat,plon,'l2m',.true.) end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET TIME STEP USING DESIRED CFL dxmin = scpx_min if(mnproc==1)print*,'TW2: dxmin',dxmin !! dx_cfl = CFL*dxmin amax = wspd(1) twave = dx_cfl/amax !! 'waves_timeres' is time [in hours] between wave readings !! - only use enough timesteps to cover this time interval nwavesteps = ceiling(waves_timeres*3600/twave) #if defined (WAVES_TESTADV) ! jump_testadv = 5 ! jump_testadv = 60 ! jump_testadv = ceiling(nwavesteps/40.0) !dint->floor jump_testadv = ceiling(20.*60/twave) !test roughly every 20mins w_test = nint(nw/2.0)! freq index to test wave advection wth_test = 4! dirn index to test wave advection ! w_test = 1 nwrite_testadv = 1 do_testadv = .false. if (mnproc==1) print*,"TWa: #t-steps, dump freq (t-steps,s)" & ,nwavesteps,jump_testadv,jump_testadv*twave #else wth_test = 2 w_test = 10 #endif/*WAVES_TESTADV flag*/ if (mnproc==1) then print*,'TWtest: twave, CFL',twave, CFL end if !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP&SCHEDULE(STATIC,jblk) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!set dave (average floe size) and WAVE MASK !! MOST THINGS ARE LOCAL SO THIS MAY NOT MAKE A DIFFERENCE !! (APART FROM ADVECTION - BUT TILING (call xctilr) !! COMMUNICATES BOUNDARY VALUES BEFORE AND AFTER THIS) do j=1,jj!j=1-nbdy,jj+nbdy do i=1,ii!i=1-nbdy,ii+nbdy if ( ICE_MASK(i,j).eq.0.0 ) then dave_store(i,j) = 0.0!!change dave_store else!!some ice present - see if need to change dave_store; ! Compute the mean floe size ! Eq. (16) Dumont et al. submitted; if ( (dmax(i,j).lt.dfloe_miz) & .and.(dmax(i,j).gt.dfloe_min) ) then ! power law distribution of Toyota et al (2011) #if defined (WAVES_FSD_RG) ! old RG method - discontinuous in Dmax ! => tau[x,y]_wav are noisy dave = dmean(fragility,xi,dfloe_min,dmax(i,j),1) #else ! smoothed version - fields are less noisy call floe_scaling_smooth(dave,fragility,xi & ,dfloe_min,dmax(i,j),1) #endif else! all floes have same length; dave = dmax(i,j) end if dave_store(i,j) = dave end if end do ! end spatial loop j end do ! end spatial loop i !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! CALC ATTENUATION COEFFICIENTS if (mnproc==1)print*, 'TW: getting attenuation coefficients' do w=1,nw period = wave_period(w) do j=1,jj do i=1,ii if (ICE_MASK(i,j).eq.1.0) then!! check ice present if (ACoption.eq.0) then alp_nd_all(i,j,w) = alpha_coeff(period,iceh(i,j)) else call attenuation_ice(alp_nd_all(i,j,w),dmpg_all(i,j,w) & ,tmp1,tmp2,period,iceh(i,j),ACoption,AGoption) !tmp1 = ag_ice, tmp2 = lam_ice if (AGoption.eq.1) then wlng_ice(i,j,w) = tmp2 if (ACoption.eq.1) then kice = 2*pi/tmp2 dmpg_all(i,j,w) = damping_RP_pert( & kice,omega(w),iceh(i,j)) ! & kice,omega(w),iceh(i,j),visc_RP,prams_RP) end if else wlng_ice(i,j,w) = wlng(w) end if !! if ((i==itest).and.(j==jtest)) then !!test damping and wavelength; print*,'' print*,'******************************************' print*,'TW_ijt: atten, damping, wavelength etc' print*,'i,j,lon,lat',i+i0,j+j0,plon(i,j),plat(i,j) print*,'period: ',2*pi/omega(w) print*,'thickness: ',iceh(i,j) print*,'damping parameter: ',visc_RP print*,'ACoption: ',ACoption print*,'AGoption: ',AGoption print*,'wavelength: ',wlng_ice(i,j,w) print*,'damping: ',dmpg_all(i,j,w) print*,'atten coeff: ',alp_nd_all(i,j,w) print*,'******************************************' print*,'' end if end if ! ag_eff = (1-icef(i,j))*wspd(w) ! & +icef(i,j)*ag_ice end if!Check for ice present end do!end i spatial loop end do!end j spatial loop end do!end w frequency loop !!******************************************************************** !!** VARIABLES INITIALISED, NOW START THE WAVE INTEGRATION LOOP: ** !!** (Time stepping) !!******************************************************************** if (mnproc==1) then print*, 'TW: Beginning time integration' print*, 'TW: current time step = ',twave,' s' print*, 'TW: Maximum CFL, amax: ', & dxmin/maxval(wspd),maxval(wspd) print*, 'TWtest: tile size: ',ii,jj end if #if defined (WAVES_TIME_INTERP) t_interp_waves = 0. #endif do n=1,nwavesteps! start time loop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !PREPARE ALL TILES FOR ADVECTION; #if defined (WAVES_TIME_INTERP) c_interp_waves = n/real(nwavesteps) t_interp_waves = t_interp_waves+twave DO_INTERP_WAVES = (t_interp_waves.ge.t_interp_period) & .or.(n.eq.nwavesteps)!do interp on last timestep if((mnproc==1).and.DO_INTERP_WAVES)then print*,' ' print*,'WAVES_TIME_INTERP' print*,'time (mins)',n*twave/60. print*,'interp freq (mins)',t_interp_period/60. print*,'t-step,interval (mins)' & ,twave/60.,twave*nwavesteps print*,' ' end if #endif #if defined (WAVES_TESTADV) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do_testadv = (mod(n,jump_testadv).eq.1) #if defined (WAVES_TIME_INTERP) do_testadv = DO_INTERP_WAVES! TESTING - overwrite "do_testadv" & "jump_testadv" ! to just test when doing interp of waves #endif do_testadv = (do_testadv.and.(nwrite_testadv.le.10)) if ((do_testadv)) then if (mnproc==1)then print*, 'waves_advect: testing advection at t-step no',n print*, 'waves_advect: time (mins)',n*twave/60. end if tmp2d = sdf_dir(:,:,wth_test,w_test) call check_adv_nc(dtime,nwrite_testadv,tmp2d) nwrite_testadv = nwrite_testadv + 1 end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif/*WAVES_TESTADV flag check*/ !!TODO !!calc integrals after advection !! - in case we want to break before advection !!mv advection to subroutine !m0_disp_2d = 0.0 !m2_disp_2d = 0.0 !m0_strain_2d = 0.0 do w=1,nw call xctilr(sdf_dir(1-nbdy,1-nbdy,1,w),1,ndir,6,6,halo_ps) !!TODO call do_advection() here !!- then integrate over frequency do wth=1,ndir sdf_ij = sdf_dir(:,:,wth,w)!thing to be advected uw = ag_eff(:,:,w)*uwave(:,:,wth) vw = ag_eff(:,:,w)*vwave(:,:,wth) !DO ADVECTION; !! previously used iceadv in mod_common_ice.F !! - now we use waveadv_weno or waveadv_lax in mod_common_wavesice.F #if defined (WAVES_LAXADV) !! Luke's LAX advection code (quicker, but needs waves of !! same speed) !! Also needs flux-limiting adv_dir = wavdir(wth) call waveadv_lax(sdf_ij,ag_eff(:,:,w),adv_dir, & scpy,scpx,twave) sdf_dir(:,:,wth,w) = sdf_ij #else !! Used to be -DWAVES_ADV_WENO flag - now default !! Luke's/Tim's adjusted advection routine (fixes land issues) call waveadv_weno(sdf_ij,uw,vw, & scuy,scux,scp2i,scp2,twave) sdf_dir(:,:,wth,w) = sdf_ij #endif/*Choice of advection scheme*/ #if defined (WAVES_TIME_INTERP) if (DO_INTERP_WAVES) then !!reset timer t_interp_waves = 0. !!overwrite waves inside overlap of current/future wave masks !!- interpolate present/future states where(WAVE_MASK2 & (1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).eq.1.0) & sdf_dir(:,:,wth,w) = sdf_out(:,:,wth,w)+ & c_interp_waves* & (sdf_dir2(:,:,wth,w)-sdf_out(:,:,wth,w)) !! testing if ((itest.gt.0).and.(jtest.gt.0)) then print*, 'TWnNc: ',dtime_waves_human, & WAVE_MASK2(itest,jtest) print*, 'TWnNc: ',n,nwavesteps,c_interp_waves print*, 'TWnNc: ',w,wth, & sdf_dir(itest,jtest,wth_test,w_test) print*, 'TWnNc: ' & ,sdf_out(itest,jtest,wth_test,w_test) & ,sdf_dir2(itest,jtest,wth_test,w_test) end if end if #endif end do! end directional loop wth end do! end spectral loop w !! ADVECTION DONE !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!do attenuation and calc integrals inside subroutine do_atten* !!do breaking inside subroutine do_breaking m0_disp_2d = 0.0 m2_disp_2d = 0.0 m0_strain_2d = 0.0 Mdir = 0.0 taux_wav_ice = 0.0 tauy_wav_ice = 0.0 taux_wav_ocn = 0.0 tauy_wav_ocn = 0.0 taux_wav = 0.0 tauy_wav = 0.0 do w=1,nw om = omega(w) !! atten_dim = 0.0 where (dave_store.gt.0.0) & & atten_dim = abs(alp_nd_all(:,:,w))*icef/dave_store !! call do_atten_simple(sdf_dir(:,:,:,w) & & ,Sfreq,taux_om,tauy_om,mth_disp_2d & & ,atten_dim,2*dmpg_all(:,:,w)*icef,ag_eff(:,:,w),twave) !! m0_disp_2d = m0_disp_2d+wt_omega(w)*Sfreq m2_disp_2d = m2_disp_2d+wt_omega(w)*Sfreq*om**2 !! tmp2d = (2*PI/wlng_ice(:,:,w))**2!!kice**2 m0_strain_2d = m0_strain_2d & & +wt_omega(w)*Sfreq*(tmp2d*iceh/2.0)**2 !! Mdir = Mdir+wt_omega(w)*mth_disp_2d !! where (ap_eff(:,:,w).gt.0.) !all the stress (momentum flux) is going into the ice !TODO add a partition in do_atten_simple !between ice and ocean taux_wav_ice = taux_wav_ice & & +wt_omega(w)*rhowtr*g*taux_om/ap_eff(:,:,w) tauy_wav_ice = tauy_wav_ice & & +wt_omega(w)*rhowtr*g*tauy_om/ap_eff(:,:,w) end where end do!w - freq loop !total ice+ocean wave stress taux_wav = taux_wav_ice+taux_wav_ocn tauy_wav = tauy_wav_ice+tauy_wav_ocn if ((itest.gt.0).and.(jtest.gt.0)) then !diagnostics for text file m0_test = m0_disp_2d(itest,jtest) m2_test = m2_disp_2d(itest,jtest) mth_test = Mdir(itest,jtest) veps_test = m0_strain_2d(itest,jtest) end if Hs = 4*sqrt(abs(m0_disp_2d)) if ((itest.gt.0).and.(jtest.gt.0)) then i = itest j = jtest print*,'TW track Hs: n,Hs',n,Hs(i,j) end if !calc mwd (Tp a dummy variable) Tp = 0.0 where (m0_disp_2d.gt.0.0) & & Tp = Mdir/m0_disp_2d Mdir = Tp !waves-from, degrees, already rotated !calc Tp Tp = 0.0 where (m2_disp_2d.gt.0.0) & & Tp = 2*PI*sqrt(abs(m0_disp_2d/m2_disp_2d)) call do_breaking(m0_strain_2d) !$OMP END PARALLEL DO end do!time loop !!******************************************************************** !!** WAVE INTEGRATION (time stepping) LOOP FINISHED: ** !!** NOW GET OUTPUTS OF SUBROUTINE AND IF DESIRED !! DUMP SOME VARIABLES INTO AN nc FILE ** !!******************************************************************** c -------------------------------------------------------------------------- ! !Diagnostic variable mwd: ! !rotate Mdir to lon-lat ref ! Mdir = -radian*(90.+Mdir)!radians,waves-to ! uwave(:,:,1) = cos(Mdir) ! vwave(:,:,1) = sin(Mdir) ! call rotate2(uwave(1-nbdy,1-nbdy,1),vwave(1-nbdy,1-nbdy,1), ! & plat,plon,'m2l',.true.) do j=1,jj do i=1,ii !rotated direction ! adv_dir = atan2(vwave(i,j,1),uwave(i,j,1))!radians,waves-to ! Mdir(i,j) = -90.-radinv*adv_dir !degrees,waves-from !finally make sure it is in [0,360) dir = Mdir(i,j) if (dir.gt.360.) then ddir = dir-360. Mdir(i,j) = dir-ceiling(ddir/360.)*360. elseif (dir.lt.0.) then ddir = -dir Mdir(i,j) = dir+ceiling(ddir/360.)*360. end if end do end do c -------------------------------------------------------------------------- !!check LAND/ICE/WAVE MASKS where (LAND_MASK.eq.1.0) MASKS = land_val where (WAVE_MASK.eq.1.0) MASKS = wave_val #if defined (WAVES_TIME_INTERP) where (WAVE_MASK2.eq.1.0) MASKS = wave_val2 #endif where (ICE_MASK.eq.1.0) MASKS = ice_val where ((dmax.gt.0).and.(dmax.lt.dfloe_miz)) & MASKS = miz_val #if defined (WAVES_CHECK_FINAL) call check_final_nc(dtime) #endif dfloe = dmax !! advect no of floes, which should be conserved !! - this is proportional to conc/Dmax**2 Nfloe = 0.0 do j=1,jj do i=1,ii !if (ICE_MASK(i,j).eq.1.0) then if (dfloe(i,j).gt.0.) then Nfloe(i,j) = icef(i,j)/dfloe(i,j)**2 end if end do end do !final outputs for ???archv_wav.* files !- will be overwritten by incident waves ! for the given output time swh_out = Hs mwp_out = Tp mwd_out = Mdir if ((itest.gt.0.).and.(jtest.gt.0)) then i = itest j = jtest txtfil = './WAVES/wavesice_out_' & //dtime_waves_human//'.txt' open (unit=txtnop,file=trim(txtfil),status='replace') !uoff+13 write(txtnop,'(i4.4,a)') i+i0 ,' # i' write(txtnop,'(i4.4,a)') j+j0 ,' # j' write(txtnop,'(e17.10,a)') plon (i,j) ,' # lon' write(txtnop,'(e17.10,a)') plat (i,j) ,' # lat' write(txtnop,'(e17.10,a)') icef (i,j) ,' # c' write(txtnop,'(e17.10,a)') iceh (i,j) ,' # h' write(txtnop,'(e17.10,a)') dfloe(i,j) ,' # Dmax' write(txtnop,'(e17.10,a)') dave_store(i,j),' # D_av' write(txtnop,'(a)') ' ' write(txtnop,'(i4.4,a)') nw ,' # nfreq' write(txtnop,'(i4.4,a)') ndir ,' # ndir' write(txtnop,'(e17.10,a)') dtheta ,' # dtheta' write(txtnop,'(a)') ' ' write(txtnop,'(e17.10,a)') Hs (i,j) ,' # Hs' write(txtnop,'(e17.10,a)') Tp (i,j) ,' # Tp' write(txtnop,'(e17.10,a)') Mdir(i,j) ,' # mwd' !!debugging diag's write(txtnop,'(a)') ' ' write(txtnop,'(e17.10,a)') mth_test ,' # mth' write(txtnop,'(e17.10,a)') m0_test ,' # m0' write(txtnop,'(e17.10,a)') m2_test ,' # m2' write(txtnop,'(e17.10,a)') veps_test ,' # veps' write(txtnop,'(e17.10,a)') young ,' # young' write(txtnop,'(e17.10,a)') radian ,' # radian' !! write(txtnop,'(a)') ' ' write(txtnop,'(e17.10,a)') taux_wav_ice(i,j) ,' # taux_wav_ice' write(txtnop,'(e17.10,a)') tauy_wav_ice(i,j) ,' # tauy_wav_ice' write(txtnop,'(e17.10,a)') taux_wav_ocn(i,j) ,' # taux_wav_ocn' write(txtnop,'(e17.10,a)') tauy_wav_ocn(i,j) ,' # tauy_wav_ocn' write(txtnop,'(e17.10,a)') taux_wav (i,j) ,' # taux_wav' write(txtnop,'(e17.10,a)') tauy_wav (i,j) ,' # tauy_wav' !wave freq's write(txtnop,'(a)') ' ' write(txtnop,'(a)') 'Omega' write(txtnop,'(a)') ' ' do w=1,nw write(txtnop,'(e15.10)') omega(w) end do write(txtnop,'(a)') ' ' !wave dir's (lon/lat) write(txtnop,'(a)') ' ' write(txtnop,'(a)') 'Wave-from dirns' write(txtnop,'(a)') ' ' do wth=1,ndir write(txtnop,'(e15.10)') wavdir(wth) end do write(txtnop,'(a)') ' ' !atten coeffs write(txtnop,'(a)') ' ' write(txtnop,'(a)') 'atten' write(txtnop,'(a)') ' ' do w=1,nw write(txtnop,'(e15.10)') & alp_nd_all(i,j,w)*icef(i,j)/dave_store(i,j) & +2*dmpg_all(i,j,w)*icef(i,j) end do write(txtnop,'(a)') ' ' !full directional spectrum write(txtnop,'(a)') ' ' write(txtnop,'(a)') 'Full spectrum' write(txtnop,'(a)') ' ' do w=1,nw do wth=1,ndir write(txtnop,'(e15.10)') sdf_dir(i,j,wth,w) end do end do close(unit=txtnop) end if call xcsync(no_flush) end subroutine waves_advect !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine check_final_nc(dtime) use mod_xc use mod_hycom_nersc, only: ncdraft!,ncerr use mod_common_wavesice use mod_common_ice,only: MIZ_MASK ! use netcdf ! use mod_year_info implicit none include "common_blocks.h" real*8,intent(in) :: dtime !! integer,parameter :: TEST_OUTPUTS = 1 !! character(len=4) :: cyy character(len=3) :: cdd,cdfrac character(len=80) :: ncfil_out !! ! integer :: ncid,dimid,var_id ! character(len=*),parameter :: tunits = ! & 'hours since 19500101T000000Z' !! & 'seconds since 1970-1-1T00:00:00Z' ! integer,parameter :: year0 = 1950 ! integer,parameter :: tfac = 24 !days to hours ! !integer,parameter :: tfac = 24*3600 !days to secs ! type(year_info) :: tt ! integer*8 :: time_val ! integer :: jday !! character(len=*),parameter :: outdir = './WAVES/out/' real, dimension(itdm,jtdm) :: fld ncfil_out = outdir//'wavesice_out_' & //dtime_waves_human//'.nc' if (mnproc==1) then print*, 'TW netcdf dump time:',dtime print*, 'TW netcdf dump time:',dtime_waves_human print*, 'TW dumping to netcdf file:',ncfil_out end if ! Dump final ice state (& c, h) in a netcdf file call xcaget(fld,plon,0) call ncdraft(trim(ncfil_out),fld,'lon','clobber') call xcaget(fld,plat,0) call ncdraft(trim(ncfil_out),fld,'lat','') ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !TODO add time dimension to nc files ! ! - this would require a new function ncdraft_time or something ! ! which would write a (idm,jdm,1) "fld" variable ! ! - would probably also require time>32767 possible ! !get time info ! call year_day(dtime,refyear,tt,yrflag) ! jday = datetojulian(tt%iyy,tt%imm,tt%idm+1,year0,1,1) ! time_val = tfac*jday+nint(tfac*(dtime-floor(dtime)))![units] to 1 Jan [year0] ! ! if (mnproc==1) then ! print*, 'TW netcdf dump time:',dtime ! print*, 'TW netcdf dump time:',dtime_waves_human ! print*, 'TW dumping to netcdf file:',ncfil_out ! print*, 'date,time:',tt%cdate,',',tt%ctime ! print*, 'time_val',time_val ! end if ! !add to nc file ! if (mnproc==1) then ! call ncerr(NF90_open(ncfil_out,NF90_WRITE,ncid)) ! !! ! call ncerr(nf90_redef(ncid)) ! call ncerr(nf90_def_dim(ncid, 'time', NF90_UNLIMITED, dimid)) ! call ncerr(NF90_DEF_VAR(ncid,'time',dimid,var_id)) ! call ncerr(NF90_PUT_ATT(ncid,var_id,'units',tunits)) ! call ncerr(NF90_PUT_ATT(ncid,var_id,'long_name', 'time')) ! call ncerr(nf90_enddef(ncid)) ! print*,'nc6' ! call ncerr(NF90_PUT_VAR(ncid,var_id,time_val))!crashes if >32767 ! !call ncerr(NF90_PUT_VAR(ncid,var_id,32767)) ! !call ncerr(NF90_PUT_VAR(ncid,var_id,32768)) ! print*,'nc7' ! !! ! call ncerr(nf90_close(ncid)) ! end if ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!ice thickness and concentration; call dump_netcdf(ncfil_out,iceh,'iceh') call dump_netcdf(ncfil_out,icef,'icef') !!dfloe is initial dmax, dmax is final state; call dump_netcdf(ncfil_out,dmax,'dmax') if (TEST_OUTPUTS) then !!dump initial Hs, Tw and mwd call dump_netcdf(ncfil_out,swh_out,'Hs0') call dump_netcdf(ncfil_out,mwp_out,'Tw0') call dump_netcdf(ncfil_out,mwd_out,'mwd0') !!initial Dmax call dump_netcdf(ncfil_out,dfloe,'dmax0') !! call xcaget(fld,MASKS,0) call ncdraft(trim(ncfil_out),fld,'MASKS','') end if #if defined (WAVES_MIZ_FRAC) !!variable for advection with dfloe miz_frac = 0.0 where((dfloe.gt.0.0).and.(dfloe.lt.dfloe_pack_init)) & miz_frac = 1.0 #endif !!dump final Hs, Tw and mwd call dump_netcdf(ncfil_out,Hs ,'Hs') call dump_netcdf(ncfil_out,Tp ,'Tw') call dump_netcdf(ncfil_out,Mdir,'mwd') !!stresses call dump_netcdf(ncfil_out,taux_wav_ice,'taux_wav_ice') call dump_netcdf(ncfil_out,tauy_wav_ice,'tauy_wav_ice') call dump_netcdf(ncfil_out,taux_wav_ocn,'taux_wav_ocn') call dump_netcdf(ncfil_out,tauy_wav_ocn,'tauy_wav_ocn') call dump_netcdf(ncfil_out,taux_wav ,'taux_wav') call dump_netcdf(ncfil_out,tauy_wav ,'tauy_wav') end subroutine check_final_nc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #if defined (WAVES_TESTADV) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine check_adv_nc(dtime,nwrite_testadv,sdf_test) use mod_xc use mod_hycom_nersc, only: ncdraft use mod_common_wavesice use mod_common_ice,only: MIZ_MASK implicit none include "common_blocks.h" real*8,intent(in) :: dtime integer,intent(in) :: nwrite_testadv real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) & ,intent(in) :: sdf_test !! integer,parameter :: TEST_OUTPUTS = 1 integer :: iyear,iday,ihour,i,j !! character(len=4) :: cyy character(len=3) :: cdd,cdfrac character(len=2) :: nc character(len=80) :: ncfil_adv !! character(len=*),parameter :: outdir2 = './WAVES/test_waveadv/' real, dimension(itdm,jtdm) :: fld ! call forday(dtime,yrflag, iyear,iday,ihour) ! write(cyy ,'(i4.4)') iyear ! write(cdd,'(i3.3)') iday-1!!start at iday=0 ! write(cdfrac,'(i3.3)') int(1.0e3*(ihour/24.0)) ! dtime_waves_human = cyy//'_'//cdd//'.'//cdfrac ncfil_adv = outdir2//'/wavesice_testadv_' & //dtime_waves_human//'.nc' if (mnproc==1) then print*, 'TW netcdf dump time:',dtime print*, 'TW netcdf dump time:',dtime_waves_human print*, 'TW dumping to netcdf file:',ncfil_adv end if if (nwrite_testadv.eq.1) then ! Dump lon,lat on 1st call only call xcaget(fld,plon,0) call ncdraft(trim(ncfil_adv),fld,'lon','clobber') call xcaget(fld,plat,0) call ncdraft(trim(ncfil_adv),fld,'lat','') end if write(nc,'(i2.2)') nwrite_testadv call dump_netcdf(ncfil_adv,sdf_test,'sdf'//nc) call dump_netcdf(ncfil_adv,Hs,'Hs'//nc) call dump_netcdf(ncfil_adv,dmax,'dmax'//nc) call dump_netcdf(ncfil_adv,Mdir,'mwd'//nc) !!stresses call dump_netcdf(ncfil_adv,taux_wav_ice,'taux_wav_ice'//nc) call dump_netcdf(ncfil_adv,tauy_wav_ice,'tauy_wav_ice'//nc) call dump_netcdf(ncfil_adv,taux_wav_ocn,'taux_wav_ocn'//nc) call dump_netcdf(ncfil_adv,taux_wav_ocn,'tauy_wav_ocn'//nc) call dump_netcdf(ncfil_adv,taux_wav ,'taux_wav' //nc) call dump_netcdf(ncfil_adv,tauy_wav ,'tauy_wav' //nc) if (mnproc==1) & print*, 'TWtest_adv: outfile = ',trim(ncfil_adv) end subroutine check_adv_nc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine do_breaking(m0_strain) use mod_xc use mod_wim_prams use mod_common_wavesice implicit none include "common_blocks.h" real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) & & ,intent(in) :: m0_strain real :: PI real :: Tcrest,wlng_crest,ccrest real :: ommin,ommax,om,dw real :: om1,lam1,lam2 real :: Pstrain,Pcrit,dave integer :: i,j,jcrest,nw logical :: BC,critter PI = 4*atan(1.0) nw = n_wavefreq ommin = minval(omega) ommax = maxval(omega) dw = (ommax-ommin)/(nw-1.) do j=1,jj do i=1,ii !!try to do breaking only if both ice and waves present critter = (ICE_MASK(i,j).eq.1.0) & & .and.(m0_strain(i,j).gt.epsc**2/2.) if (critter) then Tcrest = Tp(i,j) if (AGoption.eq.1) then!!use ice wavelength om = 2*pi/Tcrest if (om.le.ommin) then wlng_crest = wlng_ice(i,j,1) elseif (om.ge.ommax) then wlng_crest = wlng_ice(i,j,nw) else jcrest = 1+floor((om-ommin)/dw) om1 = omega(jcrest) lam1 = wlng_ice(i,j,jcrest) lam2 = wlng_ice(i,j,jcrest+1) !! ccrest = (om-om1)/dw wlng_crest = lam1+ccrest*(lam2-lam1) end if else!!water wavelength wlng_crest = gravity*Tcrest**2/(2.0*pi) end if !!update dmax and dave if breaking occurs if (wlng_crest.lt.2*dmax(i,j)) then dmax(i,j)=max(dfloe_min,0.5*wlng_crest) !!make sure dmax stays above dfloe_min !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Compute the mean floe size ! Eq. (16) Dumont et al. (2011) if ( (dmax(i,j).lt.dfloe_miz) & .and.(dmax(i,j).gt.dfloe_min) ) then ! power law distribution of Toyota et al (2011) #if defined (WAVES_FSD_RG) ! old RG method - discontinuous in Dmax ! => tau[x,y]_wav are noisy dave = dmean(fragility,xi,dfloe_min,dmax(i,j),1) #else ! smoothed version - fields are less noisy call floe_scaling_smooth(dave,fragility,xi & ,dfloe_min,dmax(i,j),1) if ((i.eq.itest).and.(j.eq.jtest)) then print*,'TW: test break',dmax(i,j),dave print*,'TW: test break (P)',fragility,xi,dfloe_min end if #endif else! all floes have same length dave = dmax(i,j) end if end if dave_store(i,j) = dave end if end do end do end subroutine do_breaking !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine do_atten_simple(Sdir & & ,Sfreq,taux_om,tauy_om,Mdir_om & & ,atten_dim,damp_dim,ag_eff,dt) use mod_xc use mod_common_wavesice, only: dfloe,dfloe_pack_init,ICE_MASK & & ,wavdir,n_wavdir implicit none integer,parameter :: ndir = n_wavdir !! real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,ndir) & ,intent(inout) :: Sdir real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) & & ,intent(out) :: Sfreq,taux_om,tauy_om,Mdir_om real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) & & ,intent(in) :: ag_eff,atten_dim,damp_dim real,intent(in) :: dt !! integer :: i,j,wth ! real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & ! & uwave,vwave real :: PI,adv_dir,S_th,tmp,alp_dim,source & ,uw,vw PI = 4.0*atan(1.0) !!weights for integral over direction Sfreq = 0.0 taux_om = 0.0 tauy_om = 0.0 Mdir_om = 0.0 do j = 1,jj do i = 1,ii if (ICE_MASK(i,j).gt.0.0) then !! atten_dim = ENERGY attenuation coeff [m^{-1}] do wth=1,ndir if (.false.) then !this gives stress on lon-lat grid adv_dir = -PI/180.0*(90.0+wavdir(wth)) uw = cos(adv_dir) vw = sin(adv_dir) else !this gives stress on x-y grid uw = uwave(i,j,wth) vw = vwave(i,j,wth) end if S_th = Sdir(i,j,wth) alp_dim = atten_dim(i,j)+damp_dim(i,j) !!stress calculation source = -alp_dim*ag_eff(i,j)*S_th!! m^{-1}*[m/s]*[m^2s] = m^2 tmp = -uw*wt_theta(wth)*source taux_om(i,j) = taux_om(i,j)+tmp tmp = -vw*wt_theta(wth)*source tauy_om(i,j) = tauy_om(i,j)+tmp !! tau_x,tau_y [m^2] need to be multiplied by !! rho_wtr*g/phase_vel [kg/m^3*m/s^2*s/m=N*s/m^4] !! and integrated over frequency as well; !! units: [m^2]*[N*s/m^4]/s = N/m^2 = Pa !! NB take '-' because here we have calc'd the stress !! ON the waves (cf Donelan et al, 2012, JGR) !! - we want the stress on the ice !!do attenuation Sdir(i,j,wth) = S_th*exp(-alp_dim*ag_eff(i,j)*dt) end do end if!!end ice check (attenuation and stress calc) !! INTEGRATE SPECTRUM OVER DIRECTION; !! NB do this in water also do wth=1,ndir Sfreq(i,j) = Sfreq(i,j)+wt_theta(wth)*abs(Sdir(i,j,wth)) Mdir_om(i,j) = Mdir_om(i,j) & +wt_theta(wth)*Sdir(i,j,wth)*wavdir(wth) end do end do!i end do!j end subroutine do_atten_simple !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function dmean(f,xi,dmin,dmax,mom) ! mean or mean square D (mom=1 or 2) ! RG method (Dumont et al, 2011) real, intent(in) :: dmin,dmax,f integer, intent(in) :: xi,mom integer :: m,mm real :: nm,dm,nsum,ndsum,r mm = 0 r = dmax/dmin do while ( r .gt. xi ) r = r/xi mm = mm+1 end do if ( mm .gt. 0 ) then nsum = 0.0 ndsum = 0.0 do m = 0,1,mm nm = (1.0-f)*(f*xi**2)**m dm = dmax/(xi**m) nsum = nsum +nm ndsum = ndsum+nm*dm**mom end do dmean = ndsum/nsum else dmean = dmin**mom end if end function dmean !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine floe_scaling_smooth(Dave,fragility,xi,dmin,dmax,mom) implicit none real, intent(out) :: Dave real, intent(in) :: dmax,fragility,dmin integer, intent(in) :: xi,mom integer :: n real :: b,Dthresh,fsd_exp,A n = mom fsd_exp = 2+log(fragility)/log(real(xi))!!power law exponent: P(d>D)=(D_min/D)^fsd_exp; b = n-fsd_exp Dthresh = 200. !! calculate from Dmax !! - uniform dist for larger floes Dave = Dmax**n !! small floes if (Dmax.le.Dmin) then Dave = Dmin**n elseif (Dmax.le.Dthresh) then !! bigger floes A = (fsd_exp*exp(fsd_exp*(log(Dmin)+log(Dmax)))) A = A/( exp(fsd_exp*log(Dmax)) & -exp(fsd_exp*log(Dmin)) ) Dave = -(A/b)*(exp(b*log(Dmin))-exp(b*log(Dmax))) end if end subroutine floe_scaling_smooth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function damping_RP_pert(k,om,h) !real function damping_RP_pert(k,om,h,visc_rp,prams_rp) use mod_wim_prams, only: rhowtr,young,poisson,visc_RP implicit none !real, dimension(3), intent(in) :: prams_RP real, intent(in) :: k,om,h!,visc_RP real :: D!,rho_wtr,young,poisson !young = prams_RP(1) !poisson = prams_RP(2) !rho_wtr = prams_RP(3) D = young*h**3/12/(1-poisson**2) damping_RP_pert = om*k**2*visc_RP/(4*D*k**5+rhowtr*om**2) end function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpha_coeff(t,h) implicit none real, intent(in) :: t,h real :: p1,p2,p3 real :: q11,q12,q13,q14,q15 real :: q21,q22,q23,q24,q25 real :: q31,q32,q33,q34,q35 ! A polynomial fit of the attenuation coefficient ! Kohout and Meylan (2008) q11 = -0.00000777 q12 = 0.00032080 q13 = -0.00437542 q14 = 0.02047559 q15 = -0.01356537 q21 = 0.00003635 q22 = -0.00153484 q23 = 0.02121709 q24 = -0.09289399 q25 = -0.03693082 q31 = -0.00004509 q32 = 0.00214484 q33 = -0.03663425 q34 = 0.26065369 q35 = -0.62474085 p1 = q11*t**4 + q12*t**3 + q13*t**2 + q14*t + q15 p2 = q21*t**4 + q22*t**3 + q23*t**2 + q24*t + q25 p3 = q31*t**4 + q32*t**3 + q33*t**2 + q34*t + q35 if ( h .lt. 0.4 ) then alpha_coeff = p1*0.4**2+p2*0.4+p3 else if ( h .gt. 5.0 ) then alpha_coeff = p1*5.0**2+p2*5.0+p3 else alpha_coeff = p1*h**2+p2*h+p3 end if if ( alpha_coeff .gt. 0.0 ) alpha_coeff = 0.0 end function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif/*if defined WAVES*/ end module mod_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) ! ! 20140721: Flags for saving waves/FSD simplified at svn Revision 330 ! WAVES_SAVE_ICEWAVES -> default (Best option) ! WAVES_SAVE_ICEWAVES2 -> WAVES_TIME_INTERP ! Default -> WAVES_NO_SAVE ! NB with WAVES_RESET_ICE flag, can choose to regularly reset the FSD, ! with a frequency of "ice_reset_freq" (set in mod_forcing_waves.F) ! ! 2015013-19: ~rev 357 ! *tidy up file: move netdcdf dumping to subroutines ! TODO: subroutines for atten and breaking ! - wrote and compiled them ! - not called or tested yet ! TODO: subroutine for advection ! *some variables (eg wavdir) moved to mod_common_wavesice.F ! to be accessible to subroutines ! TODO?: remove WAVES_TIME_INTERP option? ! - too hard to restructure around it ! - can put it back later if wanted ! TODO: look into default options outside ICE_MASK (conc<5%) ! ! 20151107: removed WAVES_NO_SAVE option ! improved WAVES_TIME_INTERP option (less freq interp) ! dump_netcdf moved to mod_common_wavesice.F ! ! 20151226: diagnostic variables mwd, wave stress (taux_wav,tauy_wav) calculated ! diagnostic text file output at (i,j)=(itest,jtest) ! -----------------------------------------------------------------------------------