module mod_forcing_waves use mod_xc implicit none logical :: DO_INIT_CHQ_MIZF contains c #if defined (WAVES) #if defined (WAVES_ONLY) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine run_waves_only!(day1,day2) use mod_common_wavesice use mod_waves_only use mod_hycom_nersc implicit none !real, intent(in) :: day1,day2 real :: day1,day2 real, parameter :: get_ice_freq = 24.0!get ice c,h daily !! "get_wave_freq" overall timestep [h] !! how often to check if its time to read waves and run the WIM !! - with WW3 frequency is every 3h (can get WW3 hourly if we want it), !! with WAM frequency is every 6h, !! so 3h (or even 1h) is a good choice real, parameter :: get_wave_freq = 3.0 ! real, parameter :: get_wave_freq = 24.0!for testing integer :: i,nsteps,ndays & ,steps_from_ice_read & ,nwav_per_day,nice_read real :: dtime0,dtime,dt_wav !!get 1st and last day from limits file: open( unit=uoff+99,file='limits') read( uoff+99,*) day1,day2 close(unit=uoff+99) day1 = abs(aint(day1)) !NB can be negative day2 = aint(day2)-1 !NB don't include last day (to be consistent with ordinary HYCOM mode) if (mnproc==1) print*,'TW (run_waves_only): day1',day1,day2 call read_ice_thick_conc(day1) call waves_set_icemask call waves_fsd_reset ndays = 1+aint(day2-day1) nwav_per_day = aint(24.0/get_wave_freq) dt_wav = 1.0/nwav_per_day nice_read = aint(get_ice_freq/get_wave_freq) nsteps = ndays*nwav_per_day steps_from_ice_read = 0 do i=1,nsteps dtime = day1+(i-1)*dt_wav dtime0 = real(floor(dtime)) !each day is separate run, so set dtime0 to start of day !Progress report call prtstatus(dtime) !!mod_hycom_nersc.F if (mnproc==1) print*, & 'TW (mod_forcing_waves.F, run_waves_only): dtime=',dtime if (steps_from_ice_read.ge.nice_read) then !! get ice conc/thickness if (mnproc==1) print*, & 'TW (run_waves_only) updating c/h: dtime=',dtime call read_ice_thick_conc(dtime) steps_from_ice_read = 0 !!reset ice mask and dfloe call waves_set_icemask call waves_fsd_reset !!reset waves sdf_out = 0.0 swh_out = 0.0 mwp_out = 0.0 end if ! check if it's time to read waves and run waves_advect call forcing_waves(dtime,dtime0) steps_from_ice_read = steps_from_ice_read+1 end do end subroutine run_waves_only !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine forcing_waves(dtime,dtime0,day2) use mod_readwaves use mod_wavesice use mod_common_ice use mod_common_wavesice use mod_hycom_nersc use netcdf implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!DECLARE VARIABLES: real*8, intent(in) :: dtime,dtime0,day2 real*8 :: dtimefrc_waves,dtimefrc2_waves,waves_rdtime & ,diff_rdtime_waves !!for determining when to read waves or to reset ice: real*8 ,save :: dtimefrc_waves_last!!this is exactly zero on first call real*8 ,save :: dtime_waves_icereset_last logical :: WAVES_READ,WAVES_RUN #if defined (WAVES_RESET_ICE) && !defined (WAVES_ONLY) real :: diff_waves_icereset,ice_reset_time real,parameter :: ice_reset_freq = 24.0!! [h] !! how often to reset dfloe - should be a multiple of wave-reading frequency) logical :: ICE_RESET #endif real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) & :: tmp_ij & ,tmp_mwd,tmp_mwp,tmp_swh real :: avg_mwd,Edir,Etot !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! waves_rdtime = waves_timeres/(24.0) dtimefrc_waves = floor(dtime/waves_rdtime)*waves_rdtime dtimefrc2_waves = dtimefrc_waves+waves_rdtime !! diff_rdtime_waves = (dtime - dtimefrc_waves_last) !! on 1st call, usual criterion may not work since !! dtimefcr_waves_last hasn't been set yet !! - so use 'DO_WAVES_INIT' from mod_common_wavesice.F if (DO_WAVES_INIT) then WAVES_READ = .true. else WAVES_READ = (diff_rdtime_waves.ge.waves_rdtime) end if WAVES_RUN = (WAVES_READ.and.(dtimefrc_waves.ne.day2)) !don't run on very last time step if (mnproc==1) then print*, 'TWtest: dtime, WAVES_READ, WAVES_RUN', & dtime, WAVES_READ,WAVES_RUN print*, 'TWtest: time for wave 1&2', & dtimefrc_waves,dtimefrc2_waves print*, 'TWtest2: time for wave 1&2', & dtimefrc_waves==dtimefrc2_waves, waves_rdtime end if if (WAVES_READ) then !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!GET INCIDENT WAVE FIELD #if defined (WAVES_TIME_INTERP) !! we have already read the waves on the 1st day !! (mod_restart.F) !! - current waves are at mwd2,mwp2,swh2 !! - save these and then get future waves if (WAVES_RUN) then tmp_mwd = mwd2 tmp_mwp = mwp2 tmp_swh = swh2 !! call read_waves(dtimefrc2_waves) !! future waves are now at mwd,mwp,swh !! - these need to be moved to mwd2,mwp2,swh2 mwd2 = mwd mwp2 = mwp swh2 = swh !!move current waves to mwd,mwp,swh mwd = tmp_mwd mwp = tmp_mwp swh = tmp_swh else !don't need to read the next bunch of waves mwd = mwd2 mwp = mwp2 swh = swh2 end if #elif ! defined (WAVES_TEST_UNIFWAVES) if (.not.DO_WAVES_INIT) then !! we have already read the waves on the 1st day !! (mod_waves_init::waves_setup() ), !! so don't need to call read_waves on 1st day call read_waves(dtimefrc_waves) end if #endif !! get wave statistics - average !! TODO: is this too slow for WAVES_DIRSPEC? !! - Don't need it really, but sometimes nice to know; tmp_ij = swh**2 call xcsum(Etot,tmp_ij,ip) tmp_ij = swh**2*mwd call xcsum(Edir,tmp_ij,ip) avg_mwd = Edir/Etot if (mnproc==1) print*, 'TWtest: avg_mwd = ',avg_mwd !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #if !defined (WAVES_ONLY) !! Define the ice mask for the WIM to use !! & set the conc/thickness for it to use !! NB This is not necessary if WAVES_ONLY is used !! - the ice mask and resetting of dfloe is done externally #if ! defined (WAVES_RESET_ICE) if (.not.DO_WAVES_INIT) then !!already done at initial time step !!- other time steps always need this call waves_set_icemask end if #else ice_reset_time = ice_reset_freq/24.0!!convert hours to days if (DO_WAVES_INIT) then ICE_RESET = .true. else call waves_set_icemask!!NB don't need this at 1st time-step diff_waves_icereset = dtimefrc_waves & -dtime_waves_icereset_last ICE_RESET = (diff_waves_icereset.ge.ice_reset_time) end if if (ICE_RESET) then !!get ice mask !!set dfloe to initial value on the ice (else 0) call waves_fsd_reset dtime_waves_icereset_last = dtimefrc_waves #if !defined (WAVES_NOSAVE) !!also reset the waves when dfloe is reset sdf_out = 0.0 swh_out = 0.0 mwp_out = 0.0 mwd_out = 0.0 #endif end if #endif/*WAVES_RESET_ICE*/ #endif/*not WAVES_ONLY*/ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Now use the incident waves: !! overwrite swh_out,mwp_out,mwd_out,sdf_out !! inside the wavemask (with swh,mwp,mwd) !! - this also sets sdf_dir2 if using -DWAVES_TIME_INTERP !! TODO add something here for -DWAVES_TEST_UNIFWAVES if(mnproc==1)print*,'TW: calling set_inc_waves' !!save original wave integrals (limited no of dirns) !!before overwriting with wave model data swh_eff = swh_out mwp_eff = mwp_out mwd_eff = mwd_out call set_inc_waves() #if ! defined (WAVES_NO_ARCHIVE) !TW: 20151019 - archiving waves before waves_advect if (mnproc.eq.1)print*, & 'TW: dumping binaries before calling'// & ' waves_advect subroutine' call archiv_waves(dtimefrc_waves,dtime0) #endif DO_WAVES_INIT = .false. dtimefrc_waves_last = dtimefrc_waves end if!!read waves check if (WAVES_RUN) then !!don't do this for the very last time step if (mnproc.eq.1)print*, & 'TW: about to go into waves_advect subroutine' call waves_advect(dtimefrc_waves) end if end subroutine forcing_waves #if defined (WAVES_MIZ_FRAC) subroutine check_miz_frac(dtime) !!routine to dump the following variables periodically: use mod_xc #if defined (WAVES) use mod_common_wavesice & , only: dfloe,dfloe_pack_init,miz_frac,miz_frac_zero #endif implicit none real, intent(in) :: dtime !!for determining when to dump ice stuff real*8 ,save :: dtime_mizf_last!!this is set to exactly zero on first call logical :: CHQ_MIZF real*8 :: mizf_timestep !!time between dumps (in days) & ,dtime_mizf & ,diff_dtime_mizf real,parameter :: mizf_freq = 1.0 mizf_timestep = mizf_freq/(24.0)!!dump frequency in days if (DO_INIT_CHQ_MIZF) then DO_INIT_CHQ_MIZF = .false. CHQ_MIZF = .false.!!don't dump at time 0 dtime_mizf_last = dtime if (mnproc==1) then print*,'TWchq_mizf: initialising check_miz_frac' end if else dtime_mizf = floor(dtime/mizf_timestep) & *mizf_timestep !! diff_dtime_mizf = (dtime_mizf - dtime_mizf_last) CHQ_MIZF = (diff_dtime_mizf.ge.mizf_timestep) end if if (CHQ_MIZF) then dtime_mizf_last = dtime if (mnproc==1) then print*,'TWchq_miz_frac: checking:',dtime!,dtime_human end if where(miz_frac.gt.miz_frac_zero) miz_frac = 0.0 dfloe = dfloe_pack_init end where else if (mnproc==1) then print*,'TWchq_miz_frac: not checking:',dtime!,dtime_human end if end if end subroutine check_miz_frac #endif subroutine waves_set_icemask use mod_common_wavesice, only: & icef,iceh,ICE_MASK use mod_common_ice, only: ficem,hicem use mod_xc implicit none !include 'common_blocks.h' !! Sets ice mask for waves !! - takes ficem,hicem from mod_common_ice.F !! - model c/h OR c/h from observations; !! - sets icef,iceh in mod_common_wavesice.F !! - c/h that the WAVES module uses; !! TODO make compatible with ICESTATE !! - take thickness as average over categories? !! - sum concentration over categories? #if defined (WAVES_ONLY) && defined (WAVES_ICEOBS) !! less coverage with SMOS than with AMSR-E, !! so use a thickness threshold to determine where the ice is !! (conservative result for MIZ location) real, parameter :: hice_min = 0.05!! [m] #endif real, parameter :: fice_min = 0.15!! should be the same as in !! m_icemodels_advect.F & mod_waves_only.F !! TODO put this into mod_common_ice.F !! but would need to be very careful ICE_MASK = 0.0 #if defined (WAVES_ONLY) && defined (WAVES_ICEOBS) !! less coverage with SMOS than with AMSR-E, !! so use a thickness threshhold where ((hicem(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).gt.hice_min) & .and.(ficem(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).gt.fice_min)) & ICE_MASK = 1.0 !! double check hice/fice if (mnproc==1) print*, & 'TW: setting ice mask for WAVES_ICEOBS' icef = ICE_MASK*ficem iceh = ICE_MASK*hicem #else !! set ICE_MASK using a concentration threshhold where (ficem(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).gt.fice_min) & ICE_MASK = 1.0 !! if (mnproc==1) print*, & 'TW: setting ice mask for WAVES in usual way' icef = ICE_MASK*ficem iceh = ICE_MASK*hicem #endif #if defined (WAVES_DBL_THICK) !! Can help to increase thickness in Fram Strait, !! as thickness is underestimated there; iceh = 2*iceh #endif end subroutine waves_set_icemask subroutine waves_fsd_reset use mod_common_ice, only:ficem use mod_common_wavesice, only: & dfloe,Nfloe,ICE_MASK,dfloe_pack_init use mod_xc implicit none if (mnproc==1) print*,'TW: in waves_fsd_reset' dfloe = 0.0 Nfloe = 0.0 where (ICE_MASK.eq.1.0) dfloe = dfloe_pack_init !!may change land value later in mod_wavesice.F where (ICE_MASK.eq.1.0) Nfloe = ficem/dfloe_pack_init**2 !!may change land value later in mod_wavesice.F end subroutine waves_fsd_reset #endif/*if defined WAVES*/ end module mod_forcing_waves ! ----------------------------------------------------------------------------------- ! 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) ! -----------------------------------------------------------------------------------