diff --git a/.github/workflows/write_logfiles.csh b/.github/workflows/write_logfiles.csh old mode 100755 new mode 100644 diff --git a/README.md b/README.md index 744753d1c..b0ac33086 100644 --- a/README.md +++ b/README.md @@ -9,10 +9,17 @@ ## The CICE Consortium sea-ice model CICE is a computationally efficient model for simulating the growth, melting, and movement of polar sea ice. Designed as one component of coupled atmosphere-ocean-land-ice global climate models, today’s CICE model is the outcome of more than two decades of community collaboration in building a sea ice model suitable for multiple uses including process studies, operational forecasting, and climate simulation. +<<<<<<< Updated upstream +======= +CICE6 includes the capabilities for measuring the impacts of waves on sea ice across a global scale. In initial testing of CICE6 it was apparent that there is no code for the propagation and attenuation of wind waves and swell through sea ice. Without waves no ice breakup will occur and all new floes are assumed to have the geometry of pancake ice, rather than using the theory of Shen et al. (2001). The Waves-in-Ice Module (WIM) module is a way to propagate waves through the ice pack in CICE6, such that preliminary runs can be completed using CICE's new FSD routine. This module is from the code used in Bennetts et al. (2017) ([https://bitbucket.org/puotila/cicewithwaves/src/master/](https://bitbucket.org/puotila/cicewithwaves/src/master/)), although some modifications have been made to make it compatible with the latest version of CICE. It was decided that CICE's ice breakup routine (Roach et al., 2018; Horvat et al., 2015) would be used within this model to keep comparisons as simple as possible. The short-coming of this decision is that wave-sea ice interactions are not coupled, we are currently investigating the effects of this. + +In our testing we found that the new floe subroutine (Shen theory) is highly influential on the floe size of the ice cover. As such, applying a north-south wave cut-off (using mean wave direction) was not sufficient. As the wave field in CICE has no memory, a large swarth of cells would frequently 'switch' between zero wave energy (waves travelling north) and non-zero wave energy (waves travelling south). This resulted in new floes being assigned to the largest floe size category in the first epoch, before being rapidly broken up by the wave field in the second epoch, causing the floe size diagnostics terms to fluctate. As such, we have applied a parameteric spreading technique paired to calculate the 2D directional spectrum, from which we intergrate over the wedge travelling southward to more accurate approximate the resulting wave energy passing into the ice cover. +>>>>>>> Stashed changes This repository contains the files and code needed to run the CICE sea ice numerical model starting with version 6. CICE is maintained by the CICE Consortium. Versions prior to v6 are found in the [CICE-svn-trunk repository](https://github.com/CICE-Consortium/CICE-svn-trunk). +<<<<<<< Updated upstream CICE consists of a top level driver and dynamical core plus the [Icepack][icepack] column physics code], which is included in CICE as a Git submodule. Because Icepack is a submodule of CICE, Icepack and CICE development are handled independently with respect to the GitHub repositories even though development and testing may be done together. [icepack]: https://github.com/CICE-Consortium/Icepack @@ -20,6 +27,19 @@ CICE consists of a top level driver and dynamical core plus the [Icepack][icepac The first point of contact with the CICE Consortium is the Consortium Community [Forum][forum]. This forum is monitored by Consortium members and also opened to the whole community. Please do not use our issue tracker for general support questions. +======= +1. The ice edge is taken to be the first equatoward cell with less than puny SIC, although this threshold parameter can be changed. +2. Read in data for these respective cells from a WW3 hindcast dataset ([https://data.csiro.au/collections/collection/CI6616v010](https://data.csiro.au/collections/collection/CI6616v010)). +3. Use a Bretscheider spectrum to define the wave spectrum, $S(\omega)$, for each cell along the edge of the wave mask (`ice_floe.F90/increment_floe_long`). +4. Apply a parametric spreading technique (cosine-sqaured) to approximate a 2D wave-energy, $E(\omega,\theta) = S(\omega)D(\theta)$ (`ice_floe.F90/increment_floe_long`). +5. Extract the component heading southward $E(\omega,\theta=\pi)$, (`ice_floe.F90/increment_floe_long`). +4. Propagate waves given the significant wave height, peak period, and mean wave direction into the ice pack. +5. Feed this data into the icepack thermodynamic routine. + +Attenuation is computed as per the observation of Meylan et al. (2014) currently, there is the code in place to implement the attenuation model presented in Williams et al. (2013a,2013b), but this has not been tested within CICE6 yet. + +For more details see https://noahday.notion.site/README-af9c753fb9ab47f2a0220bda3cfa256d. +>>>>>>> Stashed changes [forum]: https://xenforo.cgd.ucar.edu/cesm/forums/cice-consortium.146/ @@ -27,12 +47,15 @@ If you expect to make any changes to the code, we recommend that you first fork In order to incorporate your developments into the Consortium code it is imperative you follow the guidance for Pull Requests and requisite testing. Head over to our [Contributing][contributing] guide to learn more about how you can help improve CICE. +<<<<<<< Updated upstream [contributing]: https://github.com/CICE-Consortium/About-Us/wiki/Contributing ## Useful links * **CICE wiki**: https://github.com/CICE-Consortium/CICE/wiki Information about the CICE model +======= +>>>>>>> Stashed changes * **CICE Release Table**: https://github.com/CICE-Consortium/CICE/wiki/CICE-Release-Table @@ -46,5 +69,64 @@ Head over to our [Contributing][contributing] guide to learn more about how you List of resources for information about the Consortium and its repositories as well as model documentation, testing, and development. +<<<<<<< Updated upstream ## License See our [License](LICENSE.pdf) and [Distribution Policy](DistributionPolicy.pdf). +======= +Shen, H.H., Ackley, S.F., Hopkins, M.A., 2001. A conceptual model for pancake-ice formation in a wave field. Ann. Glaciol. 33, 361–367. https://doi.org/10.3189/172756401781818239 + +Williams, T.D. et al. (2013) ‘Wave–ice interactions in the marginal ice zone. Part 1: Theoretical foundations’, Ocean Modelling, 71, pp. 81–91. doi:10.1016/j.ocemod.2013.05.010. + +Williams, T.D. et al. (2013) ‘Wave–ice interactions in the marginal ice zone. Part 2: Numerical implementation and sensitivity studies along 1D transects of the ocean surface’, Ocean Modelling, 71, pp. 92–101. doi:10.1016/j.ocemod.2013.05.011. + +## ❓FAQ +- What is the domain of the module? + + The WIM is currently limited to only the Southern Hemisphere as there are more complexities with land around the Arctic. + +- What dimension are the wave spectra? + + Currently, the WIM is limited to 1D (along longitudinal lines). There is some infrastructure for moving to a 2D wave propagation scheme, although Alberto Alberello and Luke Bennetts ran into problems with the transfer of wave energy between blocks in CICE. + +- What is the wave spectrum made up of? + + Using significant wave height and a peak period a Bretschneider spectrum is defined. This spectrum is then propagated either North or South according to it's mean wave direction. + +- Where are waves propagated from? + + Wave statistics are taken from the closest ice-free cell. It uses a 1D wave scattering advection scheme, so it reads in wave direction but theses directions don't evolve over time. + +- What output does this module produce? + + Wave spectrum, peak period, significant wave height for each cell. + +- What FSD do they use? + + They use the prognostic FSD routine in CICE6 (Roach et al., 2017). However it could be possible to implement a split power law to describe the FSD following Toyota et al. 2017 (as seen in Bennetts et al., 2017) for the breakup routine. + +- Does it use wave direction? + + They do use wave direction from the closest ice free cell. The wave spectrum is 'advected' through the ice. I say 'advected' as it is really only attenuation rather than full advection. And the wave direction does not change as the waves collide with floes and other waves. + +- Relevant variable dimensions + + ```fortran + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + swh , & ! significant wave height (m) + mwd , & ! mean wave direction (Rads) + ppd ! wave peak period (s) + ``` + +- What are the attenuation settings? + + There are two kinds of attenuation present in this module, the simplest is based off of observations from Meylan et al. (2014). The more complex (floe size dependent) attenuation model has not currently been tested in CICE6 and is based off of the work by Williams et al. (2013a,2013b). + + +### Extra techinical details + +`increment_floe` is the main script in the module, it calls all other files in the `wave_ice_code` and exists within CICE code. This code first initalises: frequency min/max, angular frequencies, frequency, wave direction, initial energy spec, wavelength, wave number, dummy spectral moments, wave spectrum tolerance, spectrum passed through ice-covered ocean, sine mean wave direction for rows, values required for calculations of mean wave direction, WIM termination flags, attenuation parameters, and points in frequency and angular domains. + +This code proceeds by defining the wave mask edge and propagate waves into ice, define attenuation coefficient coefficients. Then it calculates the wavelengths and wave number by `lambda = gT^2/2pi`, `kappa = 2pi/lambda` as well as direction `th_in = -pi/2 + (lp_i-1)pi/(nth_in-1)`. Then the Bretschneider SDF function is called to initialise the wave spectrum. + +This gets the code to the point where it can update the wave spectrum and the floe sizes. This can work in both the coupled and uncoupled set ups by calling `sub_Balance` or `sub_Uncoupled` respectively. `sub_Balance` first calculates the `mass,`and `fr`, where `mass` is the scaled mass and `fr` is the flexural ridging (note that `dom` and `dth` are `d*omega` and `d(th)`). It then implements [Simpson's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule) for `dth` and calculates `lambda_ice`. Its advection scheme differs from `sub_Uncoupled` as it uses the functions `fn_IntAttn` and `fn_AvAttn` to calculate the integrated attenuation coefficient with respect to floe diameter and the averaged attenuation coefficient with respect to floe diameter respectively. This calculation takes the maximum floe size, `lambda_ice(lp_i)` and `om(lp_i)` as input. The wave attenuation is handled in the same exponential (+ cosine) scheme as in the previous subroutine (although with different coefficient values). +>>>>>>> Stashed changes diff --git a/cice.setup b/cice.setup old mode 100755 new mode 100644 diff --git a/cicecore/cicedynB/analysis/ice_history_fsd.F90 b/cicecore/cicedynB/analysis/ice_history_fsd.F90 index 43afc1e27..4f21a8737 100644 --- a/cicecore/cicedynB/analysis/ice_history_fsd.F90 +++ b/cicecore/cicedynB/analysis/ice_history_fsd.F90 @@ -163,6 +163,27 @@ subroutine init_hist_fsd_2D "floe size distribution, perimeter", & "per unit ice area", c1, c0, ns, f_fsdperim) +<<<<<<< Updated upstream +======= +! Noah Day WIM ----------------------------------------------------------------- + if (f_peak_period(1:1) /= 'x') & + call define_hist_field(n_peak_period,"peak_period","1",tstr2D, tcstr, & + "peak period of wind and swell waves", & + "from attenuated spectrum in ice", c1, c0, & + ns, f_peak_period) + if (f_mean_wave_dir(1:1) /= 'x') & + call define_hist_field(n_mean_wave_dir,"mean_wave_dir","1",tstr2D, tcstr, & + "mean wave direction of swell waves", & + "from attenuated spectrum in ice", c1, c0, & + ns, f_mean_wave_dir) + if (f_pancake_ice(1:1) /= 'x') & + call define_hist_field(n_pancake_ice,"pancake_ice","%/m",tstr2D, tcstr, & + "pancake ice concentration", & + "from aice", c1, c0, & + ns, f_pancake_ice) +! ------------------------------------------------------------------------------ + +>>>>>>> Stashed changes enddo ! nstreams diff --git a/cicecore/cicedynB/analysis/ice_history_shared.F90 b/cicecore/cicedynB/analysis/ice_history_shared.F90 index 52d268990..3389ce796 100644 --- a/cicecore/cicedynB/analysis/ice_history_shared.F90 +++ b/cicecore/cicedynB/analysis/ice_history_shared.F90 @@ -122,7 +122,7 @@ module ice_history_shared integer (kind=int_kind), parameter, public :: & nvar = 12 , & ! number of grid fields that can be written ! excluding grid vertices - nvarz = 6 ! number of category/vertical grid fields written + nvarz = 6 ! number of category/vertical grid fields written, ND: changing from 6 to 7 integer (kind=int_kind), public :: & ncat_hist , & ! number of thickness categories written <= ncat @@ -202,7 +202,12 @@ module ice_history_shared f_bounds = .true., f_NCAT = .true., & f_VGRDi = .true., f_VGRDs = .true., & f_VGRDb = .true., f_VGRDa = .true., & +<<<<<<< Updated upstream f_NFSD = .false. +======= + f_NFSD = .false., f_NFREQ = .true.!, & + !f_WFSD = .true. ! ND: adding FSD binwidth +>>>>>>> Stashed changes character (len=max_nstrm), public :: & ! f_example = 'md', & @@ -348,7 +353,12 @@ module ice_history_shared f_bounds , f_NCAT , & f_VGRDi , f_VGRDs , & f_VGRDb , f_VGRDa , & +<<<<<<< Updated upstream f_NFSD , & +======= + f_NFSD , f_NFREQ , & + !f_WFSD, & +>>>>>>> Stashed changes ! f_example , & f_hi, f_hs , & f_snowfrac, f_snowfracn, & diff --git a/cicecore/cicedynB/general/ice_floe.F90 b/cicecore/cicedynB/general/ice_floe.F90 new file mode 100644 index 000000000..7f85dcf8f --- /dev/null +++ b/cicecore/cicedynB/general/ice_floe.F90 @@ -0,0 +1,1752 @@ +!======================================================================= +! +!BOP +! +! !MODULE: ice_floe - Floe size distribution - from the cicewithwaves +! project +! +! !DESCRIPTION: +! +! !REVISION HISTORY: +! SVN:$$ +! +! authors Noah Day, using the work from +! Luke Bennetts and Siobhan O'Farrell +! date: 25-08-21 +! +! !INTERFACE: +! + module ice_floe +! +! +! !USES: +! + use icepack_kinds + use icepack_parameters, only: p01, p5, c0, c1, c2, c3, c4, c10 + use icepack_parameters, only: bignum, puny, gravit, pi + use icepack_tracers, only: nt_fsd + use icepack_warnings, only: warnstr, icepack_warnings_add, icepack_warnings_aborted + use ice_domain_size, only: ncat, nfsd, nfreq +!use ice_state ! noah day wim 004 + !use ice_constants + use ice_fileunits + use ice_read_write + use ice_restart_shared, only: lenstr, restart_dir, restart_file, & + pointer_file, runtype + use ice_communicate, only: my_task, master_task + use ice_exit, only: abort_ice + use m_prams_waveice, only: waveicedatadir, fname_ww3, WAVE_METH, ww3_lat, ww3_lon, & + ww3_dir, ww3_tm, ww3_swh, ww3_fp, ATTEN_METH, ATTEN_MODEL, attn_fac, do_coupled, & + OVERWRITE_DIRS, ww3_dir_full, ww3_swh_full, ww3_fp_full, nww3_dt + !use netcdf + use ice_constants, only: eps1, eps3 +! +!EOP +! + implicit none + + !logical (kind=log_kind) :: & + !tr_fsd, & ! if .true., use floe tracer + !restart_fsd, & ! if .true., read floe tracer restart file + !tr_wspec, & ! if .true., use wave spec tracer + !restart_spec ! if .true., read wave spec tracer restart file + +!======================================================================= + + contains + +!======================================================================= +!BOP +! +! !ROUTINE: init_floe +! +! !DESCRIPTION: +! +! Initialize ice floe tracer (call prior to reading restart data) +! +! !REVISION HISTORY: same as module +! +! !INTERFACE: +! +subroutine init_floe +! +! !USES: + use ice_arrays_column, only: floe_rad_c + use ice_state, only: trcrn + use icepack_tracers, only: nt_fsd, tr_fsd +! +!EOP +! +! integer (kind=int_kind) :: i, j + + +!real (kind=dbl_kind), dimension(:,:), intent(inout) :: & +! trcrn ! tracer array +! noah day wim 005 ------- + + call init_floe_0 + +! -------------------- + +! if (trim(runtype) == 'continue') restart_fsd = .true. + + + +! call init_floe_0 +! trcrn(i,j,nt_fsd,:,:) = c0 +! do j = 1, wavemask +! do i = 1, nx_block +! trcrn(i,j,nt_fsd,:,:) = c300 +! enddo +! enddo +! floe_rad_c(i,j,:) = c0 +! do j = 1, wavemask +! do i = 1, nx_block +! floe_rad_c(i,j,:) = c300 +! enddo +! enddo + + + end subroutine init_floe + +!======================================================================= +!BOP +! +! !ROUTINE: init_floe_0 +! +! !DESCRIPTION: +! +! Initialize ice floe tracer +! +! !REVISION HISTORY: same as module +! +! !INTERFACE: +! + subroutine init_floe_0 +! +! !USES: +! + use ice_state, only: trcrn!, nt_ifloe + use icepack_tracers, only: nt_fsd, nfsd + use ice_flux, only: ifd + use ice_blocks, only: nx_block, ny_block + use ice_constants, only: c0, c1, c2, c3, c4, c10, c30, c110, c300, c1000, p25, p01 +! +!EOP +! + integer (kind=int_kind) :: i, j +! This is broken - Noah day 20/3/22 +!trcrn(:,:,nt_fsd,:,:) = c0 +ifd(:,:,:) = c0 +do j = 1, ny_block + do i = 1, nx_block + !trcrn(:,:,nt_fsd,1,:) = c300 ! Noah Day WIM, instead of :,: at the end + ifd(i,j,:) = c300 + enddo +enddo + + +end subroutine init_floe_0 + +!======================================================================= + +!BOP +! +! !ROUTINE: increment_floe +! +! !DESCRIPTION: +! +! Increase ice floe tracer by scaled timestep length. +! +! !REVISION HISTORY: same as module +! +! !INTERFACE: +! + subroutine increment_floe (nx_block, ny_block, & + dt, tmask, Lcell, & + loc_swh, loc_ppd, loc_mwd, & + ifloe, afice, vfice, dum_wavemask, wave_spec_blk ) ! Noah Day 4/11/21 + +! !USES: +! + use m_prams_waveice + use m_waveice + use m_waveattn + use ice_fileunits, only: nu_diag + +! !INPUT/OUTPUT PARAMETERS: +! + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block ! block dimensions + +! ! time step + real (kind=dbl_kind), intent(in) :: & + dt + +! ! ice concentration & thickness (1st category) + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(in) :: & + afice, vfice + +! ! ice floe size param, significant wave height, wave peak period + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(inout) :: & + loc_swh, loc_ppd, loc_mwd, ifloe + +! ! land/boundary mask, thickness (T-cell) + logical (kind=log_kind), dimension (nx_block,ny_block), & + intent(in) :: & + tmask + + real (kind=dbl_kind), dimension (nx_block,ny_block), & + intent(in) :: & + Lcell + + integer (kind=int_kind), intent(in) :: dum_wavemask +! +! local variables +! +! ! Counters + integer (kind=int_kind) :: i, j, ij + +! ! For the wave-ice code: +! ! counter + integer :: lp_i, lp_j +! ! points in freq & angular domains + integer, parameter :: nw_in = 16, nth_in = 1 !ND: nw_in=31, nth_in=1 + !integer :: nw_in +! ! max floe size param + real(kind=dbl_kind) :: D1 +! ! cell length +! real(kind=dbl_kind), parameter :: Lcell = 25*1000d0 +! ! open ocean wave params + real(kind=dbl_kind), parameter :: fmin = 1d0/1000d0, fmax = 1d0 ! Noah Day changing from [1/16, 1/6], 1/s + ! freq min/max + real(kind=dbl_kind), parameter :: om1=2*pi*fmin, om2=2*pi*fmax ! Angular frequencies, rad/s + ! ang freqs + real(kind=dbl_kind), parameter :: om_0 = (om2 - om1)/(nw_in-1) ! rad/s + real(kind=dbl_kind), dimension(nw_in) :: om_in, T + ! freqency + real(kind=dbl_kind), dimension(nth_in) :: th_in + ! direction (not used yet!!!!) + real(kind=dbl_kind), dimension(nw_in) :: S_init_in + ! initial energy spec + + real(kind=dbl_kind), dimension(nw_in) :: lam_wtr_in, k_wtr_in + ! wavelen/num + + real(kind=dbl_kind) :: dum_sm0, dum_sm2 + ! dummy spectral moments + + real(kind=dbl_kind), parameter :: ws_tol = 1.0e-1_dbl_kind + +! ! spectrum passed through ice-covered ocean + real (kind=dbl_kind), dimension(nx_block,nw_in) :: & + wspec_row + +! ! dummy version of wspec_row to help with directionality + real (kind=dbl_kind), dimension(nx_block,nw_in) :: & + wspec_row_hld + +! ! sine mean wave directions for row + real (kind=dbl_kind), dimension(nx_block) :: & + mwd_row + +! ! array to hold values for calculation of mean wave directions + real (kind=dbl_kind), dimension(2,nx_block) :: & + mwd_hld + + + ! Noah Day 4/11/21 + real (kind=dbl_kind), dimension(nx_block,ny_block,nw_in), intent(inout) :: & + wave_spec_blk ! Noah Day 4/11/21, wave spectrum for the block + + integer, dimension(nx_block) :: tmt, tmt_hld ! WIM termination flag + +! ! attenuation + integer :: idd_alp + real(kind=dbl_kind) :: dum_alp + +!EOP +! + +!!! Define wave mask edge and propagate waves into ice: + if (cmt.ne.0) then + write(nu_diag,*) ' in increment_floe -> wavemask=', dum_wavemask + endif + +!!! Define attenuation coefficient coefficients + + if (ATTEN_MODEL.eq.1) then + write(nu_diag,*) 'opening: ', waveicedatadir , fname_alp, '...' + + open(newunit=idd_alp,file=waveicedatadir // fname_alp) + do lp_j=1,Ncheb_h+1 + do lp_i=1,Ncheb_f+1 + read(idd_alp,*) dum_alp + alp_coeffs(lp_i,lp_j)=dum_alp + end do + end do + close(idd_alp) + write(nu_diag,*) '... done' + endif + + +! Set frequency bin number equal to CICE's from ice_in + ! nw_in = nfreq + !!! Calculate wavenumbers & wavelengths + + do lp_i=1,nw_in + om_in(lp_i) = om1 + (lp_i-1)*om_0 + T(lp_i) = 2d0*pi/om_in(lp_i) + lam_wtr_in(lp_i) = gravity*(T(lp_i)**2d0)/2d0/pi + k_wtr_in(lp_i) = 2d0*pi/lam_wtr_in(lp_i) + end do + + if (nth_in.ne.1) then + do lp_i=1,nth_in + th_in(lp_i) = -pi/2 + (lp_i-1)*pi/(nth_in-1) + end do + else + th_in = 0d0 + end if + + tmt(:) = 0 ! flag to terminate routine + tmt_hld(:) = 1 ! flag to terminate routine + wspec_row(:,:) = c0 ! a dummy vector + wspec_row_hld(:,:) = c0 ! a dummy dummy vector + mwd_hld(:,:) = c0 ! another dummy dummy vector, noah day uncommented + !loc_mwd(:,:) = c0 ! another dummy dummy vector + S_init_in(:) = c0 + + !!! Begin at wavemask (only difference is initialisation) + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> starting wavemask loop' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + j = dum_wavemask + + ! loc_mwd(:,j) = pi/6d0 ! initialise, noah day uncommented + + ! A0. initialise with: Bretschneider + do i=1,nx_block + if (tmask(i,j)) then + ! initialise direction + mwd_row(i) = loc_mwd(i,j) + if (loc_swh(i,j).lt.ws_tol) then + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' NOT running wave-ice routine' + write(nu_diag,*) ' -> sig wave ht <', ws_tol + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + !write(nu_diag,*) ' -> no ice in this cell: h<', tolh + endif + if (afice(i,j).lt.tolice) then + write(nu_diag,*) ' -> no ice in this cell: c<', tolice + endif + !write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h check: swh ', loc_swh(i,j) + !write(nu_diag,*) ' -> ppd ', loc_ppd(i,j) + dum_sm0 = fn_SpecMoment(S_init_in,nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = fn_SpecMoment(S_init_in,nw_in,nth_in,om_in,th_in,2,nu_diag) + !write(nu_diag,*) ' -> swh ', 4d0*(dum_sm0**0.5d0) + !write(nu_diag,*) ' -> ppd ', & + ! 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + !write(nu_diag,*) ' -> om_in ', om_in + !write(nu_diag,*) ' -> S_init_in ', S_init_in + !write(nu_diag,*) ' -> wave_spec_blk', wave_spec_blk(i,j,:) + endif ! cmt + endif ! END COMMENT + + ! A1. Use WIM to update wave spectrum and floe sizes + + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' INPUT -> i,j =', i, j + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + write(nu_diag,*) ' -> mwd =', 180d0*mwd_row(i)/pi + endif + if (do_coupled.ne.0) then + call sub_Balance(ifloe(i,j),D1,Lcell(i,j),vfice(i,j),afice(i,j), & + nw_in,nth_in,om_in,th_in,k_wtr_in,S_init_in,wspec_row_hld(i,:),tmt(i),nu_diag) + else + call sub_Uncoupled(Lcell(i,j),vfice(i,j),afice(i,j), & + nw_in,nth_in,om_in,th_in,k_wtr_in,S_init_in,wspec_row_hld(i,:),tmt(i),nu_diag) + endif + ! Noah Day 20/3/21 ifloe(i,j) = D1 + endif ! ENDIF ws_tol + else + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> LAND' + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + endif ! ENDIF tmask(i,j) + enddo ! ENDDO i=1,nx_block + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> do wave directions', dum_wavemask + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + !print*, 'doing wave directions on wavemask loop' + + ! A2. Redistribute energy according to mean directions + ! LH Boundary: + i=1 + ! if the cell isn't land ... + if (tmask(i,j).and.tmt(i).ne.1) then + ! Calculate a dummy swh to weight directions: + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) + if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + !write(nu_diag,*) ' -> going south', 1 + !print*, ' -> going south', 1 + wspec_row(1,:) = wspec_row(1,:) + & + (1d0-sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,1) = mwd_hld(1,1) + mwd_row(1)*dum_sm0 + mwd_hld(2,1) = mwd_hld(2,1) + dum_sm0 + tmt_hld(1) = 0 + endif ! ENDIF -pi/4 going east', 1 + wspec_row(2,:) = wspec_row(2,:) + & + (sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,2) = mwd_hld(1,2) + mwd_row(1)*dum_sm0 + mwd_hld(2,2) = mwd_hld(2,2) + dum_sm0 + tmt_hld(2) = 0 + ! if wave energy needs to be advected to the west (wrap around vector) ... + elseif (tmask(nx_block,j).and.mwd_row(i).gt.pi.and.mwd_row(i).lt.3*pi/2d0) then + !write(nu_diag,*) ' -> going west', 1 + wspec_row(nx_block,:) = wspec_row(nx_block,:) + & + (sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,nx_block) = mwd_hld(1,nx_block) + mwd_row(1)*dum_sm0 + mwd_hld(2,nx_block) = mwd_hld(2,nx_block) + dum_sm0 + tmt_hld(nx_block) = 0 + endif + endif ! END IF TMASK + + ! RH Boundary: + i=nx_block + ! if the cell isn't land ... + if (tmask(i,j).and.tmt(i).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) + if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + !write(nu_diag,*) ' -> going south', nx_block + wspec_row(nx_block,:) = wspec_row(nx_block,:) + & + (1d0-sin(2d0*mwd_row(nx_block))**2d0)*wspec_row_hld(nx_block,:) + mwd_hld(1,nx_block) = mwd_hld(1,nx_block) + & + mwd_row(nx_block)*dum_sm0 + mwd_hld(2,nx_block) = mwd_hld(2,nx_block) + dum_sm0 + tmt_hld(nx_block) = 0 + endif ! ENDIF -pi/4 going south', i + wspec_row(i,:) = wspec_row(i,:) + & + (1d0-sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + mwd_hld(1,i) = mwd_hld(1,i) + mwd_row(i)*dum_sm0 + mwd_hld(2,i) = mwd_hld(2,i) + dum_sm0 + tmt_hld(i) = 0 + endif ! ENDIF -pi/4 going east', i + wspec_row(i+1,:) = wspec_row(i+1,:) + & + (sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + mwd_hld(1,i+1) = mwd_hld(1,i+1) + mwd_row(i)*dum_sm0 + mwd_hld(2,i+1) = mwd_hld(2,i+1) + dum_sm0 + tmt_hld(i+1) = 0 + elseif (tmask(i-1,j).and.mwd_row(i).gt.pi.and.mwd_row(i).lt.3*pi/2d0) then + !write(nu_diag,*) ' -> going west', i + wspec_row(i-1,:) = wspec_row(i-1,:) + & + (sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + mwd_hld(1,i-1) = mwd_hld(1,i-1) + mwd_row(i)*dum_sm0 + mwd_hld(2,i-1) = mwd_hld(2,i-1) + dum_sm0 + tmt_hld(i-1) = 0 + end if + endif ! END IF TMASK + enddo + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> update mean values', dum_wavemask + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! A3. Calculate mean parameters + do i=1,nx_block + + if (tmask(i,j)) then + if (mwd_hld(2,i).gt.c0) then + mwd_row(i) = mwd_hld(1,i)/mwd_hld(2,i) + else + mwd_row(i) = c0 + endif + tmt(i) = tmt_hld(i) + loc_mwd(i,j) = mwd_row(i) + if (tmt(i).eq.1) then + loc_swh(i,j) = 0d0 + loc_ppd(i,j) = 0d0 + else + dum_sm0 = fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,2,nu_diag) + loc_swh(i,j) = 4d0*(dum_sm0**0.5d0) + loc_ppd(i,j) = 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + wave_spec_blk(i,j,:) = wspec_row(i,:) ! Noah Day + endif + if (cmt.ne.0) then + !write(nu_diag,*) ' OUTPUT -> i,j =', i, j + !write(nu_diag,*) ' -> tmt =', tmt(i) + !write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + !write(nu_diag,*) ' -> swh =', loc_swh(i,j) + !write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + !write(nu_diag,*) ' -> mwd =', 180d0*loc_mwd(i,j)/pi + !write(nu_diag,*) ' -> wave_spec_blk(i,j,:) =', wave_spec_blk(i,j,:) + !write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + endif ! ENDIF tmask(i,j) + enddo ! ENDDO i=1,nblock + + wspec_row_hld(:,:) = c0 ! reset the dummy dummy vector + mwd_hld(:,:) = c0 ! ditto + tmt_hld(:) = 1 + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> starting wave propagation' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + do j=dum_wavemask-1,1,-1 + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> j=', j + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B1. Use WIM to update wave spectrum and floe sizes + do i=1,nx_block + !print*, 'Hi: i,j=', i, j + if (tmt(i).eq.0) then + !print*, '... tmt=0' + if (tmask(i,j)) then + if (loc_swh(i,j+1).lt.ws_tol) then + ! if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' NOT running wave-ice routine' + write(nu_diag,*) ' -> sig wave ht <', ws_tol + ! endif + tmt(i) = 1 + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + write(nu_diag,*) ' -> no ice in this cell: h<', tolh + endif + if (afice(i,j).lt.tolice) then + write(nu_diag,*) ' -> no ice in this cell: c<', tolice + endif + !write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h>>--------------------------------------------->>>' + write(nu_diag,*) ' INPUT -> i,j =', i, j + !print*, ' INPUT -> i,j =', i, j + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + !print*, ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j+1) + !print*, ' -> swh =', loc_swh(i,j+1) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j+1) + !print*, ' -> ppd =', loc_ppd(i,j+1) + write(nu_diag,*) ' -> mwd =', 180d0*mwd_row(i)/pi + !print*, ' -> mwd =', 180d0*mwd_row(i)/pi + write(nu_diag,*) ' -> afice =', afice(i,j) + ! endif ! ENDIF if (cmt.ne.0) then + if (do_coupled.ne.0) then + call sub_Balance(ifloe(i,j),D1,Lcell(i,j),vfice(i,j),afice(i,j),nw_in,& + nth_in,om_in,th_in,k_wtr_in,wspec_row(i,:),wspec_row_hld(i,:),tmt(i),nu_diag) + else + call sub_Uncoupled(Lcell(i,j),vfice(i,j),afice(i,j),nw_in,& + nth_in,om_in,th_in,k_wtr_in,wspec_row(i,:),wspec_row_hld(i,:),tmt(i),nu_diag) + endif ! ENDIF (do_coupled.ne.0) + !print*, '... done wave-ice routine' + !!ifloe(i,j) = D1 + !print*, '... set ifloe(i,j)=', ifloe(i,j), 'D1=', D1 + endif ! ENDIF ws_tol + else + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> LAND' + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + endif ! ENDIF tmask(i,j) + else + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> tmt(i)=1' + endif + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + write(nu_diag,*) ' -> no ice in this cell: h<', tolh + endif + if (afice(i,j).lt.tolice) then + write(nu_diag,*) ' -> no ice in this cell: c<', tolice + endif + !write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h do wave directions', j + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B2. Redistribute energy according to mean directions + ! LH Boundary: + i=1 + ! if the cell isn't land ... + if (tmask(1,j).and.tmt(1).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) + if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + wspec_row(1,:) = wspec_row(1,:) + & + (1d0-sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,1) = mwd_hld(1,1) + mwd_row(1)*dum_sm0 + mwd_hld(2,1) = mwd_hld(2,1) + dum_sm0 + tmt_hld(1) = 0 + endif ! ENDIF -pi/4 update mean values', j + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B3. Calculate mean parameters + do i=1,nx_block + if (tmask(i,j)) then + if (mwd_hld(2,i).gt.c0) then + mwd_row(i) = mwd_hld(1,i)/mwd_hld(2,i) + else + mwd_row(i) = c0 + endif + loc_mwd(i,j) = mwd_row(i) + tmt(i) = tmt_hld(i) + if (tmt(i).eq.1) then + loc_swh(i,j) = 0d0 + loc_ppd(i,j) = 0d0 + else + dum_sm0 = & + fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = & + fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,2,nu_diag) + loc_swh(i,j) = 4d0*(dum_sm0**0.5d0) + loc_ppd(i,j) = 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + wave_spec_blk(i,j,:) = wspec_row(i,:) ! Noah Day + endif + + ! if (cmt.ne.0) then + write(nu_diag,*) ' OUTPUT -> i,j =', i, j + write(nu_diag,*) ' -> tmt =', tmt(i) + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + write(nu_diag,*) ' -> mwd =', 180d0*loc_mwd(i,j)/pi + write(nu_diag,*) ' -> wspec_row=', wspec_row(i,:) + write(nu_diag,*) ' -> wave_spec_blk(i,j,:) =', wave_spec_blk(i,j,:) + write(nu_diag,*) '<<<---------------------------------------------<<<' + ! endif + endif ! END IF tmask(i,j) + enddo ! END do i=1,nx_block + + wspec_row_hld(:,:) = c0 ! reset the dummy dummy vector + mwd_hld(:,:) = c0 ! ditto + tmt_hld(:) = 1 + + enddo ! END do j=dum_wavemask-1,2,-1 + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> finished increment_floe' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + + + + + end subroutine increment_floe + + + !======================================================================= + + !BOP + ! + ! !ROUTINE: increment_floe_long + ! + ! !DESCRIPTION: + ! + ! Increase ice floe tracer by scaled timestep length along meridonal lines. + ! + ! !REVISION HISTORY: same as module + ! + ! !INTERFACE: + ! + subroutine increment_floe_long (nx_block, ny_block, & + dt, tmask, Lcell, & + loc_swh, loc_ppd, loc_mwd, & + ifloe, afice, vfice, dum_wavemask, dum_wavemask_vec, & + wave_spec_blk ) ! Noah Day 4/11/21 + + ! !USES: + ! + use m_prams_waveice + use m_waveice + use m_waveattn + use ice_fileunits, only: nu_diag + + ! !INPUT/OUTPUT PARAMETERS: + ! + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block ! block dimensions + + ! ! time step + real (kind=dbl_kind), intent(in) :: & + dt + + ! ! ice concentration & thickness (1st category) + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(in) :: & + afice, vfice + + ! ! ice floe size param, significant wave height, wave peak period + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(inout) :: & + loc_swh, loc_ppd, loc_mwd + +! Noah Day moving ifloe from inout to just in + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(in) :: & + ifloe ! ice floe size param + ! ! land/boundary mask, thickness (T-cell) + logical (kind=log_kind), dimension (nx_block,ny_block), & + intent(in) :: & + tmask + + real (kind=dbl_kind), dimension (nx_block,ny_block), & + intent(in) :: & + Lcell + integer (kind=int_kind), intent(in) :: dum_wavemask + integer (kind=int_kind), dimension (nx_block), intent(in) :: dum_wavemask_vec + ! + ! local variables + ! + ! ! Counters + integer (kind=int_kind) :: i, j, ij, jj + + ! Wavemask + integer (kind=int_kind) :: max_wavemask ! maximum value of the wavemask vector + + ! ! For the wave-ice code: + ! ! counter + integer :: lp_i, lp_j + ! ! points in freq & angular domains + integer, parameter :: nw_in=31, nth_in=1 + ! ! max floe size param + real(kind=dbl_kind) :: D1 + ! ! cell length + ! real(kind=dbl_kind), parameter :: Lcell = 25*1000d0 + ! ! open ocean wave params + real(kind=dbl_kind), parameter :: fmin = 1d0/1000d0, fmax = 1d0 ! ND: changing from 1/16, 1/6 + ! freq min/max + real(kind=dbl_kind), parameter :: om1=2*pi*fmin, om2=2*pi*fmax + ! ang freqs + real(kind=dbl_kind), parameter :: om_0 = (om2 - om1)/(nw_in-1) + real(kind=dbl_kind), dimension(nw_in) :: om_in, T + ! freqency + real(kind=dbl_kind), dimension(nth_in) :: th_in + ! direction (not used yet!!!!) + real(kind=dbl_kind), dimension(nw_in) :: S_init_in + ! initial energy spec + + real(kind=dbl_kind), dimension(nw_in) :: lam_wtr_in, k_wtr_in + ! wavelen/num + + real(kind=dbl_kind) :: dum_sm0, dum_sm2 + ! dummy spectral moments + + real(kind=dbl_kind), parameter :: ws_tol = 1.0e-11_dbl_kind ! ND: changing tolerance from 1.0e-1 + + ! ! spectrum passed through ice-covered ocean + real (kind=dbl_kind), dimension(nx_block,nw_in) :: & + wspec_row + + ! ! dummy version of wspec_row to help with directionality + real (kind=dbl_kind), dimension(nx_block,nw_in) :: & + wspec_row_hld + + ! ! sine mean wave directions for row + real (kind=dbl_kind), dimension(nx_block) :: & + mwd_row + + ! ! array to hold values for calculation of mean wave directions + real (kind=dbl_kind), dimension(2,nx_block) :: & + mwd_hld + + + ! Noah Day 4/11/21 + real (kind=dbl_kind), dimension(nx_block,ny_block,nw_in), intent(inout) :: & + wave_spec_blk ! Noah Day 4/11/21, wave spectrum for the block m^2s/rad + real(kind=dbl_kind) :: int_d ! Integral of the directional spectrum + + + integer, dimension(nx_block) :: tmt, tmt_hld ! WIM termination flag + + ! ! attenuation + integer :: idd_alp + real(kind=dbl_kind) :: dum_alp + + !EOP + ! + + !!! Define wave mask edge and propagate waves into ice: + if (cmt.ne.0) then + write(nu_diag,*) ' in increment_floe_long -> wavemask=', dum_wavemask_vec + endif + + !!! Define attenuation coefficient coefficients + + if (ATTEN_MODEL.eq.1) then + write(nu_diag,*) 'opening: ', waveicedatadir , fname_alp, '...' + + open(newunit=idd_alp,file=waveicedatadir // fname_alp) + do lp_j=1,Ncheb_h+1 + do lp_i=1,Ncheb_f+1 + read(idd_alp,*) dum_alp + alp_coeffs(lp_i,lp_j)=dum_alp + end do + end do + close(idd_alp) + write(nu_diag,*) '... done' + endif + + + + !!! Calculate wavenumbers & wavelengths + + do lp_i=1,nw_in + om_in(lp_i) = om1 + (lp_i-1)*om_0 + T(lp_i) = 2d0*pi/om_in(lp_i) + lam_wtr_in(lp_i) = gravity*(T(lp_i)**2d0)/2d0/pi + k_wtr_in(lp_i) = 2d0*pi/lam_wtr_in(lp_i) + end do + + if (nth_in.ne.1) then + do lp_i=1,nth_in + th_in(lp_i) = -pi/2 + (lp_i-1)*pi/(nth_in-1) + end do + else + th_in = 0d0 + end if + + tmt(:) = 0 ! flag to terminate routine + tmt_hld(:) = 1 ! flag to terminate routine + wspec_row(:,:) = c0 ! a dummy vector + wspec_row_hld(:,:) = c0 ! a dummy dummy vector + mwd_hld(:,:) = c0 ! another dummy dummy vector, noah day uncommented + !loc_mwd(:,:) = c0 ! another dummy dummy vector + S_init_in(:) = c0 ! Wave energy spectrum, m^2s/rad + + !!! Begin at wavemask (only difference is initialisation) + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> starting wavemask loop' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + + ! loc_mwd(:,j) = pi/6d0 ! initialise, noah day uncommented + + + + ! A0. initialise with: Bretschneider + do i=1,nx_block + j = dum_wavemask_vec(i) + if (j.gt.0) then + !write(nu_diag,*) 'i is:', i + !write(nu_diag,*) 'j is:', j + if (tmask(i,j)) then + ! initialise direction + mwd_row(i) = loc_mwd(i,j) + if (loc_swh(i,j).lt.ws_tol) then + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' NOT running wave-ice routine' + write(nu_diag,*) ' -> sig wave ht <', ws_tol + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + !ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + !write(nu_diag,*) ' -> no ice in this cell: h<', tolh + endif + if (afice(i,j).lt.tolice) then + !write(nu_diag,*) ' -> no ice in this cell: c<', tolice + endif + write(nu_diag,*) ' For cell (i,j)', i,j + write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h check: swh ', loc_swh(i,j) + write(nu_diag,*) ' -> ppd ', loc_ppd(i,j) + dum_sm0 = fn_SpecMoment(S_init_in,nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = fn_SpecMoment(S_init_in,nw_in,nth_in,om_in,th_in,2,nu_diag) + write(nu_diag,*) ' -> swh ', 4d0*(dum_sm0**0.5d0) + write(nu_diag,*) ' -> ppd ', & + 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + write(nu_diag,*) ' -> om_in ', om_in + write(nu_diag,*) ' -> S_init_in ', S_init_in + write(nu_diag,*) ' -> wave_spec_blk', wave_spec_blk(i,j,:) + endif ! cmt + endif ! END COMMENT + + ! A1. Use WIM to update wave spectrum and floe sizes + + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' INPUT A1 -> i,j =', i, j + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + write(nu_diag,*) ' -> mwd =', 180d0*mwd_row(i)/pi + endif + + if (do_coupled.ne.0) then + call sub_Balance(ifloe(i,j),D1,Lcell(i,j),vfice(i,j),afice(i,j), & + nw_in,nth_in,om_in,th_in,k_wtr_in,S_init_in,wspec_row_hld(i,:),tmt(i),nu_diag) + !write(nu_diag,*) '>>>--------------------------------------------->>>' + !write(nu_diag,*) 'finished sub_balance' + else + call sub_Uncoupled(Lcell(i,j),vfice(i,j),afice(i,j), & + nw_in,nth_in,om_in,th_in,k_wtr_in,S_init_in,wspec_row_hld(i,:),tmt(i),nu_diag) + endif + !ifloe(i,j) = D1 + endif ! ENDIF ws_tol + else ! Else TMASK + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> LAND' + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + endif ! ENDIF tmask(i,j) + endif ! j > 0 + enddo ! ENDDO i=1,nx_block + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> do wave directions' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + !print*, 'doing wave directions on wavemask loop' + +! ND: end of loop ------------------------------------------------------------------------------------------ + + + + + ! A2. Redistribute energy according to mean directions + ! LH Boundary: + i=1 + j = dum_wavemask_vec(i) + !write(nu_diag,*) 'j is: ', j + ! if the cell isn't land ... + if (tmask(i,j).and.tmt(i).ne.1) then + ! Calculate a dummy swh to weight directions: + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) + !if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + !write(nu_diag,*) ' -> going south', 1 + !print*, ' -> going south', 1 + wspec_row(1,:) = wspec_row(1,:) + & + (1d0-sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,1) = mwd_hld(1,1) + mwd_row(1)*dum_sm0 + mwd_hld(2,1) = mwd_hld(2,1) + dum_sm0 + tmt_hld(1) = 0 + !endif ! ENDIF -pi/4 going east', 1 + ! wspec_row(2,:) = wspec_row(2,:) + & + ! (sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + ! mwd_hld(1,2) = mwd_hld(1,2) + mwd_row(1)*dum_sm0 + ! mwd_hld(2,2) = mwd_hld(2,2) + dum_sm0 + ! tmt_hld(2) = 0 + ! ! if wave energy needs to be advected to the west (wrap around vector) ... + ! elseif (tmask(nx_block,j).and.mwd_row(i).gt.pi.and.mwd_row(i).lt.3*pi/2d0) then + ! !write(nu_diag,*) ' -> going west', 1 + ! wspec_row(nx_block,:) = wspec_row(nx_block,:) + & + ! (sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + ! mwd_hld(1,nx_block) = mwd_hld(1,nx_block) + mwd_row(1)*dum_sm0 + ! mwd_hld(2,nx_block) = mwd_hld(2,nx_block) + dum_sm0 + ! tmt_hld(nx_block) = 0 + ! endif + endif ! END IF TMASK + + ! RH Boundary: + i=nx_block + j = dum_wavemask_vec(i) + ! if the cell isn't land ... + if (tmask(i,j).and.tmt(i).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) + ! if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + !write(nu_diag,*) ' -> going south', nx_block + wspec_row(nx_block,:) = wspec_row(nx_block,:) + & + (1d0-sin(2d0*mwd_row(nx_block))**2d0)*wspec_row_hld(nx_block,:) + mwd_hld(1,nx_block) = mwd_hld(1,nx_block) + & + mwd_row(nx_block)*dum_sm0 + mwd_hld(2,nx_block) = mwd_hld(2,nx_block) + dum_sm0 + tmt_hld(nx_block) = 0 + ! endif ! ENDIF -pi/4 going south', i + wspec_row(i,:) = wspec_row(i,:) + & + (1d0-sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + mwd_hld(1,i) = mwd_hld(1,i) + mwd_row(i)*dum_sm0 + mwd_hld(2,i) = mwd_hld(2,i) + dum_sm0 + tmt_hld(i) = 0 + ! endif ! ENDIF -pi/4 going east', i + ! wspec_row(i+1,:) = wspec_row(i+1,:) + & + ! (sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + ! mwd_hld(1,i+1) = mwd_hld(1,i+1) + mwd_row(i)*dum_sm0 + ! mwd_hld(2,i+1) = mwd_hld(2,i+1) + dum_sm0 + ! tmt_hld(i+1) = 0 + ! elseif (tmask(i-1,j).and.mwd_row(i).gt.pi.and.mwd_row(i).lt.3*pi/2d0) then + ! !write(nu_diag,*) ' -> going west', i + ! wspec_row(i-1,:) = wspec_row(i-1,:) + & + ! (sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + ! mwd_hld(1,i-1) = mwd_hld(1,i-1) + mwd_row(i)*dum_sm0 + ! mwd_hld(2,i-1) = mwd_hld(2,i-1) + dum_sm0 + ! tmt_hld(i-1) = 0 + ! end if + endif ! END IF TMASK + enddo ! i = nx_block + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> update mean values', dum_wavemask_vec + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif +! ND: end of loop ------------------------------------------------------------------------------------------ + + + + + ! A3. Calculate mean parameters + do i=1,nx_block + j = dum_wavemask_vec(i) + if (tmask(i,j)) then + if (mwd_hld(2,i).gt.c0) then + mwd_row(i) = mwd_hld(1,i)/mwd_hld(2,i) + else + mwd_row(i) = c0 + endif + tmt(i) = tmt_hld(i) + loc_mwd(i,j) = mwd_row(i) + if (tmt(i).eq.1) then + loc_swh(i,j) = 0d0 + loc_ppd(i,j) = 0d0 + else + dum_sm0 = fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,2,nu_diag) + loc_swh(i,j) = 4d0*(dum_sm0**0.5d0) + loc_ppd(i,j) = 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + wave_spec_blk(i,j,:) = wspec_row(i,:) ! Noah Day + endif + if (cmt.ne.0) then + write(nu_diag,*) ' OUTPUT A3 -> i,j =', i, j + write(nu_diag,*) ' -> tmt =', tmt(i) + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + write(nu_diag,*) ' -> mwd =', 180d0*loc_mwd(i,j)/pi + write(nu_diag,*) ' -> wave_spec_blk(i,j,:) =', wave_spec_blk(i,j,:) + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + endif ! ENDIF tmask(i,j) + enddo ! ENDDO i=1,nblock + +! ND: end of loop ------------------------------------------------------------------------------------------ + + + + + wspec_row_hld(:,:) = c0 ! reset the dummy dummy vector + mwd_hld(:,:) = c0 ! ditto + tmt_hld(:) = 1 + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> starting wave propagation' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) 'MINIMUM: ', minval(dum_wavemask_vec) + write(nu_diag,*) 'max dum_wavemask is:', maxval(dum_wavemask_vec) + endif + +! if (maxval(dum_wavemask_vec).gt.0) then +! ! There is ice +! max_wavemask = maxval(dum_wavemask_vec) +! else +! ! No ice +! max_wavemask = 50 +! end if ! maxval + +max_wavemask = dum_wavemask +!write(nu_diag,*) 'max dum_wavemask is:', max_wavemask + do jj=1,max_wavemask!maxval(dum_wavemask_vec)!dum_wavemask-1,2,minval(dum_wavemask_vec)+1!-1 + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> jj=', jj + write(nu_diag,*) 'dum_wavemask is:', dum_wavemask + + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B1. Use WIM to update wave spectrum and floe sizes + do i=1,nx_block + j = dum_wavemask_vec(i) - jj!(dum_wavemask_vec(i) - jj) ! wavemask - (wavemask-j) = j + !write(nu_diag,*) '>>>--------------------------------------------->>>' + !write(nu_diag,*) 'B1. j is:', j + !write(nu_diag,*) ' i is:', i + !write(nu_diag,*) 'dum_wavemask_vec is:', dum_wavemask_vec(i) + !write(nu_diag,*) 'jj is:', jj + !write(nu_diag,*) '<<<---------------------------------------------<<<' + if (j.gt.0) then + !print*, 'Hi: i,j=', i, j + if (tmt(i).eq.0) then + !print*, '... tmt=0' + if (tmask(i,j)) then + if (loc_swh(i,j+1).lt.ws_tol) then + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' NOT running wave-ice routine' + write(nu_diag,*) ' -> sig wave ht <', ws_tol + endif + tmt(i) = 1 + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + !ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + write(nu_diag,*) ' -> no ice in this cell: h<', tolh + write(nu_diag,*) ' -> coords', i,j + endif + if (afice(i,j).lt.tolice) then + write(nu_diag,*) ' -> no ice in this cell: c<', tolice + write(nu_diag,*) ' -> coords', i,j + endif + !write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h>>--------------------------------------------->>>' + write(nu_diag,*) ' INPUT B1 -> i,j =', i, j + !print*, ' INPUT -> i,j =', i, j + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + !print*, ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j+1) + !print*, ' -> swh =', loc_swh(i,j+1) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j+1) + !print*, ' -> ppd =', loc_ppd(i,j+1) + write(nu_diag,*) ' -> mwd =', 180d0*mwd_row(i)/pi + !print*, ' -> mwd =', 180d0*mwd_row(i)/pi + endif ! ENDIF if (cmt.ne.0) then + if (do_coupled.ne.0) then + call sub_Balance(ifloe(i,j),D1,Lcell(i,j),vfice(i,j),afice(i,j),nw_in,& + nth_in,om_in,th_in,k_wtr_in,wspec_row(i,:),wspec_row_hld(i,:),tmt(i),nu_diag) + else + call sub_Uncoupled(Lcell(i,j),vfice(i,j),afice(i,j),nw_in,& + nth_in,om_in,th_in,k_wtr_in,wspec_row(i,:),wspec_row_hld(i,:),tmt(i),nu_diag) + endif ! ENDIF (do_coupled.ne.0) + !print*, '... done wave-ice routine' + ! Noah Day 20/3/22 ifloe(i,j) = D1 + !print*, '... set ifloe(i,j)=', ifloe(i,j), 'D1=', D1 + endif ! ENDIF ws_tol + else ! tmask + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> LAND' + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + tmt(i) = 1 + endif ! ENDIF tmask(i,j) + else ! tmt + if (cmt.ne.0) then + write(nu_diag,*) '>>>--------------------------------------------->>>' + write(nu_diag,*) ' i, j = ', i, j + write(nu_diag,*) ' -> tmt(i)=1' + endif + if (afice(i,j).lt.tolice.or.vfice(i,j).lt.tolh) then + !ifloe(i,j) = floe_sz_pancake + if (cmt.ne.0) then + if (vfice(i,j).lt.tolh) then + write(nu_diag,*) ' -> no ice in this cell: h<', tolh + endif + if (afice(i,j).lt.tolice) then + write(nu_diag,*) ' -> no ice in this cell: c<', tolice + endif + write(nu_diag,*) ' -> ifloe(i,j)=', ifloe(i,j) + endif ! END IF COMMENT + endif ! END IF h0 + enddo ! ENDDO i=1,nblock + + wspec_row(:,:) = c0 ! reset dummy vector + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> do wave directions 2'!, j + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B2. Redistribute energy according to mean directions + ! LH Boundary: + i=1 + j = dum_wavemask_vec(i) - jj!dum_wavemask_vec(i) - (dum_wavemask_vec(i) - jj) ! wavemask - (wavemask-j) = j + !write(nu_diag,*) ' Redistribute energy according to mean directions ' + !write(nu_diag,*) ' LH Boundary: j: ', j + ! if the cell isn't land ... + if (j.gt.0) then + if (tmask(1,j).and.tmt(1).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) +! if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + ! South + wspec_row(1,:) = wspec_row(1,:) + & + (1d0-sin(2d0*mwd_row(1))**2d0)*wspec_row_hld(1,:) + mwd_hld(1,1) = mwd_hld(1,1) + mwd_row(1)*dum_sm0 + mwd_hld(2,1) = mwd_hld(2,1) + dum_sm0 + tmt_hld(1) = 0 +! endif ! ENDIF -pi/4 0 + ! RH Boundary: + i=nx_block + j = dum_wavemask_vec(i) - jj!dum_wavemask_vec(i) - (dum_wavemask_vec(i) - jj) ! wavemask - (wavemask-j) = j + ! if the cell isn't land ... + if (j.gt.0) then + if (tmask(nx_block,j).and.tmt(nx_block).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) +! if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + ! South + wspec_row(nx_block,:) = wspec_row(nx_block,:) + & + (1d0-sin(2d0*mwd_row(nx_block))**2d0)*wspec_row_hld(nx_block,:) + mwd_hld(1,nx_block) = mwd_hld(1,nx_block) + & + mwd_row(nx_block)*dum_sm0 + mwd_hld(2,nx_block) = mwd_hld(2,nx_block) + dum_sm0 + tmt_hld(nx_block) = 0 +! endif ! ENDIF -pi/4 0 + + ! loop the inner cells... + do i=2,nx_block-1 + j = dum_wavemask_vec(i) - jj!dum_wavemask_vec(i) - (dum_wavemask_vec(i) - jj) ! wavemask - (wavemask-j) = j + if (j.gt.0) then + if (tmask(i,j).and.tmt(i).ne.1) then + dum_sm0 = & + fn_SpecMoment(wspec_row_hld(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm0 = 4d0*(dum_sm0**0.5d0) +! if (mwd_row(i).gt.3d0*pi/4d0.and.mwd_row(i).lt.5d0*pi/4d0) then + ! South + wspec_row(i,:) = wspec_row(i,:) + & + (1d0-sin(2d0*mwd_row(i))**2d0)*wspec_row_hld(i,:) + mwd_hld(1,i) = mwd_hld(1,i) + mwd_row(i)*dum_sm0 + mwd_hld(2,i) = mwd_hld(2,i) + dum_sm0 + tmt_hld(i) = 0 +! endif ! ENDIF -pi/4 0 + end do ! ENDDO i=2,nx_block-1 + + if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> update mean values'!, j + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooo' + endif + + ! B3. Calculate mean parameters + do i=1,nx_block + j = dum_wavemask_vec(i) - jj!dum_wavemask_vec(i) - (dum_wavemask_vec(i) - jj) ! wavemask - (wavemask-j) = j + if (j.gt.0) then + if (tmask(i,j)) then + if (mwd_hld(2,i).gt.c0) then + mwd_row(i) = mwd_hld(1,i)/mwd_hld(2,i) + else + mwd_row(i) = c0 + endif + loc_mwd(i,j) = mwd_row(i) + tmt(i) = tmt_hld(i) + if (tmt(i).eq.1) then + loc_swh(i,j) = 0d0 + loc_ppd(i,j) = 0d0 + else + dum_sm0 = & + fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,0,nu_diag) + dum_sm2 = & + fn_SpecMoment(wspec_row(i,:),nw_in,nth_in,om_in,th_in,2,nu_diag) + loc_swh(i,j) = 4d0*(dum_sm0**0.5d0) + loc_ppd(i,j) = 2d0*pi*((dum_sm0/dum_sm2)**0.5d0) + wave_spec_blk(i,j,:) = wspec_row(i,:) ! Noah Day + endif + + if (cmt.ne.0) then + write(nu_diag,*) ' OUTPUT B3 -> i,j =', i, j + write(nu_diag,*) ' -> tmt =', tmt(i) + write(nu_diag,*) ' -> ifloe =', ifloe(i,j) + write(nu_diag,*) ' -> swh =', loc_swh(i,j) + write(nu_diag,*) ' -> ppd =', loc_ppd(i,j) + write(nu_diag,*) ' -> mwd =', 180d0*loc_mwd(i,j)/pi + write(nu_diag,*) ' -> wspec_row=', wspec_row(i,:) + write(nu_diag,*) ' -> wave_spec_blk(i,j,:) =', wave_spec_blk(i,j,:) + write(nu_diag,*) '<<<---------------------------------------------<<<' + endif + endif ! END IF tmask(i,j) + endif ! j > 0 + enddo ! END do i=1,nx_block + + wspec_row_hld(:,:) = c0 ! reset the dummy dummy vector + mwd_hld(:,:) = c0 ! ditto + tmt_hld(:) = 1 + +enddo ! END do j = dum_wavemask-1,2,-1 +! ND: end of loop ------------------------------------------------------------------------------------------ + + + + + +if (cmt.ne.0) then + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' + write(nu_diag,*) ' -> finished increment_floe_long' + write(nu_diag,*) 'oooooooooooooooooooooooooooooooooooooooooooooooo' +endif + + + + + + end subroutine increment_floe_long + +!======================================================================= +!---! these subroutines write/read Fortran unformatted data files .. +!======================================================================= +! +!BOP +! +! !IROUTINE: write_restart_floe - dumps all fields required for restart +! +! !INTERFACE: +! + subroutine write_restart_floe(filename_spec) +! +! !DESCRIPTION: +! +! Dumps all values needed for restarting +! +! !REVISION HISTORY: +! +! author Elizabeth C. Hunke, LANL +! +! !USES: +! + use ice_domain_size + use ice_calendar, only: msec, mmonth, mday, myear, istep1, & + timesecs, idate, year_init ! Noah Day WIM, removed , time_forc + use ice_state +! +! !INPUT/OUTPUT PARAMETERS: +! + character(len=char_len_long), intent(in), optional :: filename_spec + +!EOP +! + integer (kind=int_kind) :: & + i, j, k, n, it, iblk, & ! counting indices + iyear, imonth, iday ! year, month, day + + character(len=char_len_long) :: filename + + logical (kind=log_kind) :: diag + + ! construct path/file + if (present(filename_spec)) then + filename = trim(filename_spec) + else + iyear = myear + year_init - 1 + imonth = mmonth + iday = mday + + write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & + restart_dir(1:lenstr(restart_dir)), & + restart_file(1:lenstr(restart_file)),'.floe.', & + iyear,'-',mmonth,'-',mday,'-',msec + end if + + ! begin writing restart data + !call ice_open(nu_dump_floe,filename,0) + + if (my_task == master_task) then + !write(nu_dump_floe) istep1,timesecs!,time_forc + write(nu_diag,*) 'Writing ',filename(1:lenstr(filename)) + endif + + diag = .true. + + !----------------------------------------------------------------- + + do n = 1, ncat + !call ice_write(nu_dump_floe,0,trcrn(:,:,nt_fsd,n,:),'ruf8',diag) + enddo + + !if (my_task == master_task) close(nu_dump_floe) + + end subroutine write_restart_floe + +!======================================================================= +!BOP +! +! !IROUTINE: read_restart_floe - reads all fields required for restart +! +! !INTERFACE: +! + subroutine read_restart_floe(filename_spec) +! +! !DESCRIPTION: +! +! Reads all values needed for an ice floe restart +! +! !REVISION HISTORY: +! +! author S O'Farrell +! +! !USES: +! + use ice_domain_size + use ice_calendar, only: msec, mmonth, mday, myear, istep1, & + timesecs, idate, year_init ! Noah Day WIM, removed , time_forc + use ice_state +! +! !INPUT/OUTPUT PARAMETERS: +! + character(len=char_len_long), intent(in), optional :: filename_spec + +!EOP +! + integer (kind=int_kind) :: & + i, j, k, n, it, iblk, & ! counting indices + iyear, imonth, iday ! year, month, day + + character(len=char_len_long) :: & + filename, filename0, string1, string2 + + logical (kind=log_kind) :: & + diag + + if (my_task == master_task) then + open(nu_rst_pointer,file=pointer_file) + read(nu_rst_pointer,'(a)') filename0 + filename = trim(filename0) + close(nu_rst_pointer) + + ! reconstruct path/file + n = index(filename0,trim(restart_file)) + if (n == 0) call abort_ice('ifloe restart: filename discrepancy') + string1 = trim(filename0(1:n-1)) + string2 = trim(filename0(n+lenstr(restart_file):lenstr(filename0))) + write(filename,'(a,a,a,a)') & + string1(1:lenstr(string1)), & + restart_file(1:lenstr(restart_file)),'.floe.', & + string2(1:lenstr(string2)) + endif ! master_task + + !call ice_open(nu_restart_floe,filename,0) + + if (my_task == master_task) then + !read(nu_restart_floe) istep1,timesecs!,time_forc + write(nu_diag,*) 'Reading ',filename(1:lenstr(filename)) + endif + + diag = .true. + + !----------------------------------------------------------------- + + do n = 1, ncat + !call ice_read(nu_restart_floe,0,trcrn(:,:,nt_fsd,n,:),'ruf8',diag) + enddo + + !if (my_task == master_task) close(nu_restart_floe) + + end subroutine read_restart_floe + + +end module ice_floe + +!======================================================================= diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index bcc7305ff..6a1abf510 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -16,7 +16,7 @@ module ice_flux use ice_fileunits, only: nu_diag use ice_blocks, only: nx_block, ny_block use ice_domain_size, only: max_blocks, ncat, max_nstrm, nilyr - use ice_constants, only: c0, c1, c5, c10, c20, c180 + use ice_constants, only: c0, c1, c5, c10, c20, c100, c180 use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_tracer_flags, icepack_query_tracer_indices diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index a71e6dd17..720119211 100755 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -33,8 +33,13 @@ module ice_forcing use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, & timer_bound, timer_forcing use ice_arrays_column, only: oceanmixed_ice, restore_bgc +<<<<<<< Updated upstream use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c15, c20, & c180, c360, c365, c1000, c3600 +======= + use ice_constants, only: c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c12, c15, c20, & + c180, c360, c365, c1000, c3600, c100 +>>>>>>> Stashed changes use ice_constants, only: p001, p01, p1, p2, p25, p5, p6 use ice_constants, only: cm_to_m use ice_constants, only: field_loc_center, field_type_scalar, & @@ -50,7 +55,14 @@ module ice_forcing get_forcing_atmo, get_forcing_ocn, get_wave_spec, & read_clim_data, read_clim_data_nc, & interpolate_data, interp_coeff_monthly, & +<<<<<<< Updated upstream read_data_nc_point, interp_coeff +======= + read_data_nc_point, interp_coeff, init_wave_spec, & + init_wave_spec_usr, init_wave_spec_init, check, init_wave_spec_long, & + init_wave_block + ! Noah Day WIM, adding init_wave_spec, init_wave_spec_usr, init_wave_spec_init, check +>>>>>>> Stashed changes integer (kind=int_kind), public :: & ycycle , & ! number of years in forcing cycle, set by namelist @@ -135,8 +147,14 @@ module ice_forcing wave_spec_file,& ! file name for wave spectrum oceanmixed_file ! file name for ocean forcing data +<<<<<<< Updated upstream integer (kind=int_kind), parameter :: & nfld = 8 ! number of fields to search for in forcing file +======= + integer (kind=int_kind), parameter :: & + nfld = 8, & ! number of fields to search for in forcing file + ndfld = 4 ! ND number of ocean input variables +>>>>>>> Stashed changes ! as in the dummy atm (latm) real (kind=dbl_kind), parameter, public :: & @@ -146,7 +164,8 @@ module ice_forcing frcidf = 0.17_dbl_kind ! frac of incoming sw in near IR diffuse band real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: & - ocn_frc_m ! ocn data for 12 months + ocn_frc_m, & ! ocn data for 12 months + ocn_frc_m_access ! ocn data for access logical (kind=log_kind), public :: & restore_ocn ! restore sst if true @@ -170,7 +189,7 @@ module ice_forcing ! PRIVATE: real (dbl_kind), parameter :: & - mixed_layer_depth_default = c20 ! default mixed layer depth in m + mixed_layer_depth_default = c100 ! default mixed layer depth in m logical (kind=log_kind), parameter :: & local_debug = .false. ! local debug flag @@ -212,6 +231,7 @@ subroutine alloc_forcing topmelt_data(nx_block,ny_block,2,max_blocks,ncat), & botmelt_data(nx_block,ny_block,2,max_blocks,ncat), & ocn_frc_m(nx_block,ny_block, max_blocks,nfld,12), & ! ocn data for 12 months + ocn_frc_m_access(nx_block,ny_block, max_blocks,ndfld,12), & ! ND: ocn data for 12 months from ACCESS topmelt_file(ncat), & botmelt_file(ncat), & stat=ierr) @@ -375,7 +395,15 @@ subroutine init_forcing_ocn(dt) if (trim(ocn_data_type) == 'clim') then - sss_file = trim(ocn_data_dir)//'/sss.mm.100x116.da' ! gx3 only + if (nx_global == 320) then ! gx1 + sss_file = trim(ocn_data_dir)//'/sss.mm.100x116.da' ! gx3 only + elseif (nx_global == 360) then ! access-om2 + sss_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) + elseif (nx_global == 1440) then ! access-om2-025 + sss_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) + else ! gx3 + sss_file = trim(ocn_data_dir)//'/sss.mm.100x116.da' ! gx3 only + endif if (my_task == master_task) then write (nu_diag,*) ' ' @@ -424,6 +452,10 @@ subroutine init_forcing_ocn(dt) if (nx_global == 320) then ! gx1 sst_file = trim(ocn_data_dir)//'/sst_clim_hurrell.dat' + elseif (nx_global == 360) then ! access-om2-10 + sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) + elseif (nx_global == 1440) then ! access-om2-025 + sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) else ! gx3 sst_file = trim(ocn_data_dir)//'/sst.mm.100x116.da' endif @@ -498,6 +530,10 @@ subroutine init_forcing_ocn(dt) call ocn_data_hycom_init endif + if (trim(ocn_data_type) == 'access-om2') then + call ocn_data_access_init + endif + end subroutine init_forcing_ocn !======================================================================= @@ -709,6 +745,8 @@ subroutine get_forcing_ocn (dt) call ocn_data_oned elseif (trim(ocn_data_type) == 'hycom') then ! call ocn_data_hycom(dt) + elseif (trim(ocn_data_type) == 'access-om2') then + call ocn_data_access(dt) !MHRI: NOT IMPLEMENTED YET endif @@ -1489,6 +1527,10 @@ subroutine file_year (data_file, yr) i = index(data_file,'.nc') - 5 tmpname = data_file write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' + elseif (trim(ocn_data_type)=='access-om2') then ! netcdf + i = index(data_file,'.nc') - 5 + tmpname = data_file + write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' else ! LANL/NCAR naming convention i = index(data_file,'.dat') - 5 tmpname = data_file @@ -2457,6 +2499,10 @@ subroutine JRA55_data character(len=64) :: fieldname !netcdf field name character (char_len_long) :: uwind_file_old character(len=*), parameter :: subname = '(JRA55_data)' +<<<<<<< Updated upstream +======= + !write(nu_diag,*) 'local_debug: ', local_debug +>>>>>>> Stashed changes if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start' @@ -2491,7 +2537,8 @@ subroutine JRA55_data !------------------------------------------------------------------- uwind_file_old = uwind_file - if (uwind_file /= uwind_file_old .and. my_task == master_task) then + !if (uwind_file /= uwind_file_old .and. my_task == master_task) then + if (my_task == master_task) then write(nu_diag,*) subname,' reading forcing file = ',trim(uwind_file) endif @@ -2607,6 +2654,24 @@ subroutine JRA55_data write(nu_diag,*) subname,'fdbg c12intp = ',c1intp,c2intp endif + + !vmin = global_minval(Tair_data(:,:,n1,:),distrb_info,tmask) + ! vmax = global_maxval(Tair_data(:,:,n1,:),distrb_info,tmask) + !if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg Tair_data',vmin,vmax + + ! ND: NOT NEEDED + !do iblk = 1, nblocks + ! Remove negative Tair_data + !do j = 1, ny_block + ! do i = 1, nx_block + ! if (Tair_data(i,j,n1,iblk).lt.c0) Tair_data(i,j,n1,iblk) = 185.15 ! ND: If Temp (K) < 0 set it to last value + ! enddo + !enddo + !enddo + !vmin = global_minval(Tair_data(:,:,n1,:),distrb_info,tmask) + !vmax = global_maxval(Tair_data(:,:,n1,:),distrb_info,tmask) + !if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg Tair_data',vmin,vmax + ! Interpolate call interpolate_data (Tair_data, Tair) call interpolate_data (uatm_data, uatm) @@ -2626,6 +2691,9 @@ subroutine JRA55_data do j = 1, ny_block do i = 1, nx_block if (aice(i,j,iblk) > p1) Tair(i,j,iblk) = min(Tair(i,j,iblk), Tffresh+p1) + !if (Tair(i,j,iblk).lt.185.15) Tair(i,j,iblk) = 185.15 ! ND: 185.15 is coldest ever recorded + !if (j.gt.(ny_block/2)) uatm_data(i,j,n1,iblk) = c0 + !if (j.gt.(ny_block/2)) vatm_data(i,j,n1,iblk) = c0 enddo enddo @@ -5402,3 +5470,1351 @@ end module ice_forcing !======================================================================= +<<<<<<< Updated upstream +======= +! +! !INPUT/OUTPUT PARAMETERS: +! +! +!EOP +! + integer (kind=int_kind) :: lp,lp_b,lp_i,lp_j + integer (kind=int_kind), intent(in) :: dum_wavemask + integer, intent(in) :: N_lon !, DUM_METH + real(kind=dbl_kind), dimension(N_lon), intent(in) :: dum_swh, dum_fp, dum_mwd + + integer, dimension(1) :: dumlonloc, ind_lon_ww3 !, dumlatloc +! real(kind=dbl_kind) :: dumlat, dumlon, dumswh + integer :: ind_lon + + integer (kind=int_kind), intent(in) :: & + iblk ! block index + + ! local variables + + integer (kind=int_kind) :: & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + i, j, len ! horizontal indices + + + type (block) :: & + this_block ! block information for current block + + do lp_b=1,nblocks + do lp_i=1,nx_block + do lp_j=1,ny_block + swh(lp_i,lp_j,lp_b) = c0 + ppd(lp_i,lp_j,lp_b) = c0 + mwd(lp_i,lp_j,lp_b) = c0 + enddo + enddo + enddo + + !!!!!!!!!!!!!!!!!!!!! WW3 !!!!!!!!!!!!!!!!!!!!! +! ind_lon = floor(mod(c180*TLON(1,dum_wavemask,1)/pi,c360)) + + ! Find the location at which the the ww3 grid is = 0. +!write(nu_diag,*) ' ww3_lon', SHAPE(TRANSPOSE(ww3_lon)) +!write(nu_diag,*) ' ww3_lon', SIZE(ww3_lon) +!write(nu_diag,*) ' abs(ww3_lon)', SIZE(abs(ww3_lon)) +! ind_lon_ww3 =minval(ww3_lon, mask=ww3_lon.gt.0) +! write(nu_diag,*) ' ind_lon_ww3', ind_lon_ww3 +! ind_lon_ww3 =minval(abs(ww3_lon)) +! write(nu_diag,*) ' ind_lon_ww3', ind_lon_ww3 +! write(nu_diag,*) 'minloc(abs(ww3_lon-c0),dim=2)', minloc(abs(ww3_lon-c0),dim=2) +! write(nu_diag,*) 'minloc(abs(ww3_lon-c0),dim=1)', minloc(abs(ww3_lon-c0),dim=1) + + +! NOAH DAY START +ind_lon_ww3 = minloc(abs(ww3_lon-c0),dim=1) ! find the difference in grids +! CICE starts at lon = 0, ww3 does not. + +this_block = get_block(blocks_ice(iblk),iblk) +ilo = this_block%ilo +ihi = this_block%ihi +jlo = this_block%jlo +jhi = this_block%jhi + + + !do lp_b=1,nblocks + +len = size(dum_swh) +if (cmt.ne.0) then + write(nu_diag,*) ' <--------------- block info: -------------------->' + write(nu_diag,*) ' nblocks: ', nblocks + write(nu_diag,*) ' ind_lon: ', ind_lon + write(nu_diag,*) ' nx_block: ', nx_block + !write(nu_diag,*) ' dum_swh: ', dum_swh + !write(nu_diag,*) ' swh: ', SHAPE(swh(:,dum_wavemask,:)) + write(nu_diag,*) ' N_lon: ', N_lon + write(nu_diag,*) ' dum_swh: ', SHAPE(dum_swh) + write(nu_diag,*) ' ind_lon_ww3', SHAPE(ind_lon_ww3) + write(nu_diag,*) '<-------------------------------------------------->' +end if +ind_lon = ind_lon_ww3(1)-1 ! 36th element is 359.93750000000000 + +i = ilo ! index for nx_block +j = 1 ! block index + +do lp=1,N_lon + if (lp+ind_lon.gt.N_lon) then ! if index exceeds the length of the data + if (dum_swh(lp+ind_lon-N_lon).gt.puny.and.dum_fp(lp+ind_lon-N_lon).gt.puny) then + swh(i,dum_wavemask,j) = dum_swh(lp+ind_lon-N_lon) + ppd(i,dum_wavemask,j) = c1/dum_fp(lp+ind_lon-N_lon) + !if (pi*dum_mwd(lp+ind_lon-N_lon)/c180.gt.pi/2.and.pi*dum_mwd(lp+ind_lon-N_lon)/c180.lt.3*pi/2) then + mwd(i,dum_wavemask,j) = pi*dum_mwd(lp+ind_lon-N_lon)/c180 + ! else + ! mwd(i,dum_wavemask,j) = c0 + ! end if + endif + else ! feed the horizontally translated data into CICE + if (dum_swh(lp+ind_lon).gt.puny.and.dum_fp(lp+ind_lon).gt.puny) then + swh(i,dum_wavemask,j) = dum_swh(lp+ind_lon) + ppd(i,dum_wavemask,j) = c1/dum_fp(lp+ind_lon) + !if (pi*dum_mwd(lp+ind_lon)/c180.gt.2*pi/3.and.pi*dum_mwd(lp+ind_lon)/c180.lt.4*pi/3) then + mwd(i,dum_wavemask,j) = pi*dum_mwd(lp+ind_lon)/c180 + !else + ! mwd(i,dum_wavemask,j) = c0 + !end if + endif + end if ! lp+ind_lon + i = i + 1 ! tick up until i = ihi + if (mod(i,ihi+1).eq.0) then + i = ilo ! reset nx_block index, 1 is a ghost cell + j = j + 1 ! set for the next block + end if +end do + +call ice_HaloUpdate (swh, halo_info, & + field_loc_center, field_type_scalar) +call ice_HaloUpdate (ppd, halo_info, & + field_loc_center, field_type_scalar) +call ice_HaloUpdate (mwd, halo_info, & + field_loc_center, field_type_scalar) + + +! GHOST CELLS +! Replicate the value for ilo into the ghost cells +!do j=2,nblocks +! this_block = get_block(blocks_ice(j),j) +! ilo = this_block%ilo +! ihi = this_block%ihi +! jlo = this_block%jlo +! jhi = this_block%jhi +! do i = 1,ilo-1 +! swh(i,dum_wavemask,j) = swh(ilo,dum_wavemask,j) +! ppd(i,dum_wavemask,j) = ppd(ilo,dum_wavemask,j) +! mwd(i,dum_wavemask,j) = mwd(ilo,dum_wavemask,j) +! end do +!end do + + +! ! now do it for the block below +!! j = j-1 ! block below +!! this_block = get_block(blocks_ice(j),j) +!! ilo = this_block%ilo +!! ihi = this_block%ihi +!! jlo = this_block%jlo +!! jhi = this_block%jhi +!! do i = ihi+1,nx_block +!! swh(i,dum_wavemask,j) = swh(ilo,dum_wavemask,j+1) ! now j+1 = j +!! ppd(i,dum_wavemask,j) = ppd(ilo,dum_wavemask,j+1) +!! mwd(i,dum_wavemask,j) = mwd(ilo,dum_wavemask,j+1) +!! end do +!!end do + + +! at end +!do j=1,nblocks +! this_block = get_block(blocks_ice(j),j) +! ilo = this_block%ilo +! ihi = this_block%ihi +! jlo = this_block%jlo +! jhi = this_block%jhi +! do i = ihi,nx_block +! swh(i,dum_wavemask,j) = swh(ihi,dum_wavemask,j) +! ppd(i,dum_wavemask,j) = ppd(ihi,dum_wavemask,j) +! mwd(i,dum_wavemask,j) = mwd(ihi,dum_wavemask,j) +! end do +!end do + + + + + + +!if (iblk.lt.12) then +!do lp_b=1,nblocks + +!do lp=1,nx_block!-ind_lon ! loop from 1 to dum_wavemask position +! if (dum_swh(lp+90).gt.puny.and.dum_fp(lp).gt.puny) then +! swh(lp,40,lp_b) = c1 + + !write(nu_diag,*) 'SWH in initialisation:', dum_swh(lp+ind_lon) + !write(nu_diag,*) 'puny is: ', puny + !if (dum_swh(lp).gt.puny.and.dum_fp(lp).gt.puny) then + !if (lp_b.eq.1) swh(:,dum_wavemask,lp_b) = c1 + !if (lp_b.eq.2) swh(nx_block+1-lp,dum_wavemask,lp_b) = c2 + !if (lp_b.eq.3) swh(nx_block+1-lp,dum_wavemask,lp_b) = c3 + !if (lp_b.eq.4) swh(nx_block+1-lp,dum_wavemask,lp_b) = c4 + !if (lp_b.eq.5) swh(nx_block+1-lp,dum_wavemask,lp_b) = c5 + !if (lp_b.eq.6) swh(nx_block+1-lp,dum_wavemask,lp_b) = c6 + !if (lp_b.eq.7) swh(nx_block+1-lp,dum_wavemask,lp_b) = c7 + !if (lp_b.eq.8) swh(nx_block+1-lp,dum_wavemask,lp_b) = c8 + !if (lp_b.eq.9) swh(nx_block+1-lp,dum_wavemask,lp_b) = c9 + !if (lp_b.eq.10) then + ! swh(nx_block+1-lp,dum_wavemask,lp_b) = c10 + !else + ! swh(lp,dum_wavemask,lp_b) = c1 + ! endif + + + !swh(nx_block+1-lp,dum_wavemask,lp_b) = !dum_swh(lp)!+ind_lon) +! ppd(lp,dum_wavemask,lp_b) = c1/dum_fp(lp+ind_lon) +!write(nu_diag,*) ' swh :', dum_swh(lp), dum_swh(lp+ind_lon) + + + !if (dum_mwd(lp+ind_lon).gt.c180) then + ! mwd(lp,dum_wavemask,lp_b) = pi*(dum_mwd(lp+ind_lon)-c360)/c180 + !else +! mwd(lp,dum_wavemask,lp_b) = pi*dum_mwd(lp+ind_lon)/c180 + !endif + ! else + ! swh(lp,dum_wavemask,lp_b) = c0 + !ppd(lp,dum_wavemask,lp_b) = c0 + !mwd(lp,dum_wavemask,lp_b) = c0 +! end if ! lp gt +!end do +!end do +!end if ! iblk + + + ! enddo ! end lp=1,360-ind_lon ! above the wave mask + ! do lp=363-ind_lon,nx_block + ! if (dum_swh(lp-360+ind_lon).gt.puny.and.dum_fp(lp-360+ind_lon).gt.puny) then + ! swh(lp,dum_wavemask,lp_b) = c1!dum_swh(lp-360+ind_lon) + ! ppd(lp,dum_wavemask,lp_b) = c1/dum_fp(lp-360+ind_lon) + !if (dum_mwd(lp-360+ind_lon).gt.c180) then + ! mwd(lp,dum_wavemask,lp_b) = pi*(dum_mwd(lp+ind_lon)-c360)/c180 + !else + ! mwd(lp,dum_wavemask,lp_b) = pi*dum_mwd(lp-360+ind_lon)/c180 + !endif + ! else + ! swh(lp,dum_wavemask,lp_b) = c0 + !ppd(lp,dum_wavemask,lp_b) = c0 + ! mwd(lp,dum_wavemask,lp_b) = c0 + !endif + ! enddo ! end lp=361-ind_lon,nx_block + ! enddo ! end lp_b=1,nblocks + + ! NOAH DAY END + + +!WIM START + !write(nu_diag,*) ' ind_lon : ', floor(mod(c180*TLON(dum_wavemask,1,1)/pi,c360)) + 90 + !write(nu_diag,*) ' ind_lon : ', floor(mod(c180*TLON(1,dum_wavemask,1)/pi,c360)) + 90 + +! do lp_b=1,nblocks +! do lp=1,360-ind_lon + !write(nu_diag,*) 'SWH in initialisation:', dum_swh(lp+ind_lon) + !write(nu_diag,*) 'puny is: ', puny +! if (dum_swh(lp+ind_lon).gt.puny.and.dum_fp(lp+ind_lon).gt.puny) then +! swh(lp,dum_wavemask,lp_b) = dum_swh(lp+ind_lon) +! ppd(lp,dum_wavemask,lp_b) = c1/dum_fp(lp+ind_lon) +!write(nu_diag,*) ' swh :', dum_swh(lp), dum_swh(lp+ind_lon) + + + !if (dum_mwd(lp+ind_lon).gt.c180) then + ! mwd(lp,dum_wavemask,lp_b) = pi*(dum_mwd(lp+ind_lon)-c360)/c180 + !else +! mwd(lp,dum_wavemask,lp_b) = pi*dum_mwd(lp+ind_lon)/c180 + !endif + !else + !swh(lp,dum_wavemask,lp_b) = c0 + !ppd(lp,dum_wavemask,lp_b) = c0 + !mwd(lp,dum_wavemask,lp_b) = c0 +! endif +! enddo ! end lp=1,360-ind_lon +! do lp=363-ind_lon,nx_block +! if (dum_swh(lp-360+ind_lon).gt.puny.and.dum_fp(lp-360+ind_lon).gt.puny) then +! swh(lp,dum_wavemask,lp_b) = dum_swh(lp-360+ind_lon) +! ppd(lp,dum_wavemask,lp_b) = c1/dum_fp(lp-360+ind_lon) + !if (dum_mwd(lp-360+ind_lon).gt.c180) then + ! mwd(lp,dum_wavemask,lp_b) = pi*(dum_mwd(lp+ind_lon)-c360)/c180 + !else +! mwd(lp,dum_wavemask,lp_b) = pi*dum_mwd(lp-360+ind_lon)/c180 + !endif + !else + !swh(lp,dum_wavemask,lp_b) = c0 + !ppd(lp,dum_wavemask,lp_b) = c0 + !mwd(lp,dum_wavemask,lp_b) = c0 +! endif +! enddo ! end lp=361-ind_lon,nx_block +! enddo ! end lp_b=1,nblocks +! WIM END +!len = 1 +!do lp_b=1,nblocks + +!end do + + end subroutine init_wave_spec + +!======================================================================= +!BOP +! +! !ROUTINE: init_wave_spec_usr +! +! !DESCRIPTION: +! +! Initialize wave spectrum: significant wave height and peak period (call prior to reading restart data) +! +! !REVISION HISTORY: +! +! author: L Bennetts, U Adelaide +! +! !INTERFACE: +! + subroutine init_wave_spec_usr(dum_wavemask) +! +! !USES: +! + use ice_domain, only: nblocks + use ice_flux, only: swh, ppd, mwd + use m_prams_waveice, only: pi + use ice_grid, only: tlon, tlat +! +! !INPUT/OUTPUT PARAMETERS: +! +! +!EOP +! + integer (kind=int_kind) :: lp, lp_b,lp_i,lp_j + integer (kind=int_kind), intent(in) :: dum_wavemask + + do lp_b=1,nblocks + do lp_i=1,nx_block + do lp_j=1,ny_block + swh(lp_i,lp_j,lp_b) = c0 + ppd(lp_i,lp_j,lp_b) = c0 + mwd(lp_i,lp_j,lp_b) = c0 + enddo + enddo + enddo + + do lp_b=1,nblocks + do lp=1,nx_block + swh(lp,dum_wavemask,lp_b) = c3 ! 0.28_dbl_kind ! + ppd(lp,dum_wavemask,lp_b) = c10 + mwd(lp,dum_wavemask,lp_b) = pi!7d0*pi/6d0 + enddo + enddo + + end subroutine init_wave_spec_usr + !======================================================================= + !BOP + ! + ! !ROUTINE: init_wave_spec_long + ! + ! !DESCRIPTION: + ! + ! Initialize wave spectrum per longitude: significant wave height and peak period (call prior to reading restart data) + ! + ! !REVISION HISTORY: + ! + ! author: N Day, U Adelaide + ! + ! !INTERFACE: + ! + subroutine init_wave_spec_long(dum_wavemask_vec,dum_swh,dum_fp,dum_mwd,N_lon,N_lat,iblk) + ! + ! !USES: + ! + use ice_flux, only: swh, ppd, mwd + use m_prams_waveice, only: pi, ww3_lat, ww3_lon, ww3_dir, cmt + use netcdf + use ice_grid, only: tlon, tlat + use icepack_parameters, only: puny ! Noah Day WIM + use ice_blocks, only: block, get_block, nx_block, ny_block + use ice_domain, only: blocks_ice, nblocks, halo_info ! Noah Day WIM + use ice_boundary, only: ice_HaloUpdate ! Noah Day WIM + use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global ! Noah Day WIM + use ice_communicate, only: my_task, master_task + + ! + ! !INPUT/OUTPUT PARAMETERS: + ! + ! + !EOP + ! + integer (kind=int_kind) :: lp,lp_b,lp_i,lp_j + integer (kind=int_kind), dimension(nx_block), intent(in) :: dum_wavemask_vec + !integer (kind=int_kind) :: dum_wavemask + integer, intent(in) :: N_lon, N_lat !, DUM_METH + real(kind=dbl_kind), dimension(N_lon,N_lat), intent(in) :: dum_swh, dum_fp, dum_mwd + + real(kind=dbl_kind) :: dumlon + + integer, dimension(1) :: dumlonloc, ind_lon_ww3, ind_lon_dim !, dumlatloc + ! real(kind=dbl_kind) :: dumlat, dumlon, dumswh + integer :: ind_lon + + integer (kind=int_kind), intent(in) :: & + iblk ! block index + + ! local variables + + integer (kind=int_kind) :: & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + i, j, len ! horizontal indices + + + type (block) :: & + this_block ! block information for current block + + do lp_b=1,nblocks + do lp_i=1,nx_block + do lp_j=1,ny_block + swh(lp_i,lp_j,lp_b) = c0 + ppd(lp_i,lp_j,lp_b) = c0 + mwd(lp_i,lp_j,lp_b) = c0 + enddo + enddo + enddo + + !!!!!!!!!!!!!!!!!!!!! WW3 !!!!!!!!!!!!!!!!!!!!! + + ! NOAH DAY START + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + if (cmt.ne.0) then + write(nu_diag,*) ' c180*TLON(1,dum_wavemask_vec(1),1)/pi:', c180*TLON(1,dum_wavemask_vec(1),1)/pi + write(nu_diag,*) ' mod(c180*TLON(1,dum_wavemask_vec(1),1)/pi,c360):', mod(c180*TLON(1,dum_wavemask_vec(1),1)/pi,c360) + write(nu_diag,*) ' floor(mod(c180*TLON(1,dum_wavemask_vec(1),1)/pi,c360)) :', floor(mod(c180*TLON(1,dum_wavemask_vec(1),1)/pi,c360)) + write(nu_diag,*) ' minloc(abs(ww3_lon-ind_lon),dim=2):', minloc(abs(ww3_lon-ind_lon),dim=2) + write(nu_diag,*) ' minloc(abs(ww3_lon-ind_lon),dim=1):', minloc(abs(ww3_lon-ind_lon),dim=1) + !write(nu_diag,*) ' ind_lon:', ind_lon + !write(nu_diag,*) ' ind_lon_ww3:', ind_lon_ww3 + !write(nu_diag,*) ' ind_lon_ww3(1):', ind_lon_ww3(1) + write(nu_diag,*) ' ww3_lon :', ww3_lon + write(nu_diag,*) ' <--------------- block info: -------------------->' + write(nu_diag,*) ' nblocks: ', nblocks + !write(nu_diag,*) ' ind_lon: ', ind_lon + !write(nu_diag,*) ' ind_lon_dim: ', ind_lon_dim + write(nu_diag,*) ' nx_block: ', nx_block + write(nu_diag,*) ' ilo: ', ilo + write(nu_diag,*) ' ihi: ', ihi + write(nu_diag,*) ' dum_wavemask_vec(1): ', dum_wavemask_vec(1) + write(nu_diag,*) ' block number: ', iblk + write(nu_diag,*) ' nx_global: ', nx_global + write(nu_diag,*) ' dumlon: ', dumlon + write(nu_diag,*) ' dumlonloc: ', dumlonloc + write(nu_diag,*) ' dumlonloc(1): ', dumlonloc(1) + write(nu_diag,*) ' tlon loc: ', c180*tlon(dumlonloc,dum_wavemask_vec(1),iblk)/pi + write(nu_diag,*) ' ww3 lon loc: ', c180*ww3_lon(dumlonloc,dum_wavemask_vec(1))/pi + write(nu_diag,*) ' my task: ', my_task + write(nu_diag,*) ' master_task: ', master_task + write(nu_diag,*) ' TLON bound', mod(c180*TLON(ilo,dum_wavemask_vec(1),iblk)/pi,c360), mod(c180*TLON(ihi,dum_wavemask_vec(1),iblk)/pi,c360) + write(nu_diag,*) ' ww3 bound (1,180)', mod(c180*ww3_lon(1,dum_wavemask_vec(1))/pi,c360), mod(c180*ww3_lon(180,dum_wavemask_vec(1))/pi,c360) + write(nu_diag,*) ' TLON shape', SHAPE(tlon) + write(nu_diag,*) ' ww3_lon shape', SHAPE(ww3_lon) + !write(nu_diag,*) ' dum_swh: ', dum_swh + write(nu_diag,*) ' swh: ', SHAPE(swh) + write(nu_diag,*) ' N_lon: ', N_lon + write(nu_diag,*) ' dum_swh: ', SHAPE(dum_swh) + !write(nu_diag,*) ' ind_lon_ww3', SHAPE(ind_lon_ww3) + write(nu_diag,*) '<-------------------------------------------------->' + end if + +! Find the minimum block longitude +dumlon = mod(c180*TLON(1,dum_wavemask_vec(1),iblk)/pi +c360 ,c360) +! Find the closest match on the WW3 grid +dumlonloc=minloc(abs(mod(ww3_lon+c360,c360)-dumlon),dim=1) +!write(nu_diag,*) 'BLOCK NUMBER : ', iblk +!write(nu_diag,*) 'We are seeking to find the longitude: ', dumlon +!write(nu_diag,*) 'dumlonloc: ', dumlonloc(1) +!write(nu_diag,*) 'Error, abs(ww3_lon-dumlon): ', abs(mod(ww3_lon+c360,c360)-dumlon) + + do lp_b = 1,nblocks ! should be iblk + do lp_i=ilo,ihi + ! Adapting the WW3 index for the current block/processor + i = lp_i-1 + dumlonloc(1) + !write(nu_diag,*) ' lp_i: ', lp_i + !write(nu_diag,*) ' i: ', i + + if (i.lt.N_lon) then + if (dum_wavemask_vec(lp_i).gt.0) then ! If there is a valid wave mask then read in waves + j = dum_wavemask_vec(lp_i) + !write(nu_diag,*) 'ww3_lon: ',mod(c180*ww3_lon(i,j)/pi + c360,c360) + !write(nu_diag,*) 'TLON: ', mod(c180*TLON(lp_i,j,iblk)/pi + c360,c360) + + if (dum_swh(i,j).gt.puny.and.dum_fp(i,j).gt.puny) then + !write(nu_diag,*) ' CICE lon: ', mod(c180*TLON(lp_i,dum_wavemask_vec(1),iblk)/pi + c360,c360) + !write(nu_diag,*) ' WW3 lon: ', mod(ww3_lon(i,1) + 360,c360) + swh(lp_i,j,lp_b) = dum_swh(i,j) + ppd(lp_i,j,lp_b) = c1/dum_fp(i,j) ! Converting to peak period + mwd(lp_i,j,lp_b) = dum_mwd(i,j)*c2*pi/c360 ! Converting to Radians + else + swh(lp_i,j,lp_b) = c0 + ppd(lp_i,j,lp_b) = c0 + mwd(lp_i,j,lp_b) = c0 + endif + endif + else + ! Exceeded dimension, time to start over from 1. + if (dum_wavemask_vec(lp_i).gt.0) then + j = dum_wavemask_vec(lp_i) + if (dum_swh(i-N_lon,j).gt.puny.and.dum_fp(i-N_lon,j).gt.puny) then + + !write(nu_diag,*) ' CICE lon: ', mod(c180*TLON(lp_i,dum_wavemask_vec(1),iblk)/pi + c360,c360) + !write(nu_diag,*) ' WW3 lon: ', mod(ww3_lon(i-N_lon,1) + 360,c360) + swh(lp_i,j,lp_b) = dum_swh(i-N_lon,j) + ppd(lp_i,j,lp_b) = c1/dum_fp(i-N_lon,j) ! Converting to peak period + mwd(lp_i,j,lp_b) = dum_mwd(i-N_lon,j)*c2*pi/c360 ! Converting to Radians + else + swh(lp_i,j,lp_b) = c0 + ppd(lp_i,j,lp_b) = c0 + mwd(lp_i,j,lp_b) = c0 + endif + endif + end if + end do ! wave mask +enddo ! block + + + + call ice_HaloUpdate (swh, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloUpdate (ppd, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloUpdate (mwd, halo_info, & + field_loc_center, field_type_scalar) + + end subroutine init_wave_spec_long + + + !======================================================================= + !BOP + ! + ! !ROUTINE: init_wave_block + ! + ! !DESCRIPTION: + ! + ! Distribute the read in wave statistics to each block + ! + ! !REVISION HISTORY: + ! + ! author: N Day, U Adelaide + ! + ! !INTERFACE: + ! + subroutine init_wave_block(dum_swh,dum_fp,dum_mwd,N_lon,N_lat) + ! + ! !USES: + ! + use ice_flux, only: swh, ppd, mwd + use m_prams_waveice, only: pi, ww3_lat, ww3_lon, ww3_dir, cmt, & + ww3_swh_blk, ww3_fp_blk, ww3_dir_blk + use netcdf + use ice_grid, only: tlon, tlat + use icepack_parameters, only: puny ! Noah Day WIM + use ice_blocks, only: block, get_block, nx_block, ny_block + use ice_domain, only: blocks_ice, nblocks, halo_info ! Noah Day WIM + use ice_boundary, only: ice_HaloUpdate ! Noah Day WIM + use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global ! Noah Day WIM + + ! + ! !INPUT/OUTPUT PARAMETERS: + ! + ! + !EOP + ! + integer (kind=int_kind) :: lp,lp_b,lp_i,lp_j + !integer (kind=int_kind) :: dum_wavemask + integer, intent(in) :: N_lon, N_lat + real(kind=dbl_kind), dimension(N_lon,N_lat), intent(in) :: dum_swh, dum_fp, dum_mwd + + integer, dimension(1) :: dumlonloc, ind_lon_ww3 !, dumlatloc + ! real(kind=dbl_kind) :: dumlat, dumlon, dumswh + integer :: ind_lon + + ! local variables + + integer (kind=int_kind) :: & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + i, j, len ! horizontal indices + + + type (block) :: & + this_block ! block information for current block + + ! Allocate the block variables + allocate(ww3_swh_blk(nx_block,ny_block,nblocks)) + allocate(ww3_fp_blk(nx_block,ny_block,nblocks)) + allocate(ww3_dir_blk(nx_block,ny_block,nblocks)) + + do lp_b=1,nblocks + do lp_i=1,nx_block + do lp_j=1,ny_block + ww3_swh_blk(lp_i,lp_j,lp_b) = c0 + ww3_fp_blk(lp_i,lp_j,lp_b) = c0 + ww3_dir_blk(lp_i,lp_j,lp_b) = c0 + enddo + enddo + enddo + + !!!!!!!!!!!!!!!!!!!!! WW3 !!!!!!!!!!!!!!!!!!!!! + + ! Initialise the indices for the full data + i = 1 + j = 1 + do lp_b = 1,nblocks + + this_block = get_block(blocks_ice(lp_b),lp_b) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + !if (cmt.ne.0) then + write(nu_diag,*) ' <--------------- block info: -------------------->' + write(nu_diag,*) ' nblocks: ', nblocks + write(nu_diag,*) ' ind_lon: ', ind_lon + write(nu_diag,*) ' nx_block: ', nx_block + write(nu_diag,*) ' ilo: ', ilo + write(nu_diag,*) ' ihi: ', ihi + !write(nu_diag,*) ' dum_swh: ', dum_swh + write(nu_diag,*) ' swh: '!, SHAPE(swh(:,dum_wavemask,:)) + write(nu_diag,*) ' dum_swh: ', SHAPE(dum_swh) + write(nu_diag,*) ' ind_lon_ww3', SHAPE(ind_lon_ww3) + write(nu_diag,*) '<-------------------------------------------------->' + !end if + + ! Loop over the blocks to distribute the map + do lp_i = ilo,ihi + do lp_j = jlo,jhi + i = lp_i - 1 ! Start at 1 + j = lp_j + ww3_swh_blk(lp_i,lp_j,lp_b) = dum_swh(i,j) + ww3_fp_blk(lp_i,lp_j,lp_b) = dum_fp(i,j) + ww3_dir_blk(lp_i,lp_j,lp_b) = dum_mwd(i,j) + enddo + enddo + + + + enddo ! lp_b + + call ice_HaloUpdate (ww3_swh_blk, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloUpdate (ww3_fp_blk, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloUpdate (ww3_dir_blk, halo_info, & + field_loc_center, field_type_scalar) + + end subroutine init_wave_block + + +!======================================================================= +!BOP +! +! !ROUTINE: check +! +! !DESCRIPTION: +! +! For WW3 netcdf data +! +! !REVISION HISTORY: +! +! author: L Bennetts, U Adelaide +! +! !INTERFACE: +! + subroutine check(status) + use netcdf + integer, intent ( in) :: status + if(status /= nf90_noerr) then + write(nu_diag,*) trim(nf90_strerror(status)) + stop "Stopped" + end if + end subroutine check + + + + + +! Noah Day WIM ======================================================================= + + subroutine ocn_data_access_init + +! Reads ocean output for ACCESS-OM2 +! +! List of ocean forcing fields: Note that order is important! +! (order is determined by field list in vname). +! +! For ocean mixed layer-----------------------------units +! +! 1 sst------temperature---------------------------(K) +! 2 sss------salinity------------------------------(psu=ppt) +! 3 u--------surface u current---------------------(m/s) +! 4 v--------surface v current---------------------(m/s) +! +! Fields 3, 4 are on the U-grid; 1, 2 are +! on the T-grid. + +! authors: Noah Day, UAdl (adapted from ocn_data_ncar_init) + + use ice_blocks, only: nx_block, ny_block + use ice_domain_size, only: max_blocks + use ice_domain, only: nblocks, distrb_info + use ice_flux, only: sss, sst, Tf, uocn, vocn, frzmlt + use ice_grid, only: hm, tmask, umask + use ice_global_reductions, only: global_minval, global_maxval + use ice_restart_shared, only: runtype + use ice_calendar, only: myear + use icepack_parameters, only: puny + use ice_read_write, only: ice_read_nc_uv +#ifdef USE_NETCDF + use netcdf +#endif + + integer (kind=int_kind) :: & + n , & ! field index + m , & ! month index + nrec, & ! record number for direct access + nbits + + character(char_len) :: & + vname(ndfld) ! variable names to search for in file + data vname / & + 'sst', 'sss', 'u', 'v' / + + integer (kind=int_kind) :: & + ! fid , & ! file id + dimid ! dimension id + + integer (kind=int_kind) :: & + status , & ! status flag + nlat , & ! number of longitudes of data + nlon , & ! number of latitudes of data + nzlev ! z level of currents + + integer (kind=int_kind) :: & + i, j, iblk , & ! horizontal indices + fid ! file id for netCDF file + + character (char_len) :: & + fieldname ! field name in netcdf file + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + work1 + + real (kind=dbl_kind) :: & + vmin, vmax + + character(len=*), parameter :: subname = '(ocn_data_access_init)' + + if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start' + + if (my_task == master_task) then + + write (nu_diag,*) 'WARNING: evp_prep calculates surface tilt' + write (nu_diag,*) 'WARNING: stress from geostrophic currents,' + write (nu_diag,*) 'WARNING: not data from ocean forcing file.' + write (nu_diag,*) 'WARNING: Alter ice_dyn_evp.F90 if desired.' + + if (restore_ocn) write (nu_diag,*) & + 'SST restoring timescale = ',trestore,' days' + + sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) ! not just sst + + !--------------------------------------------------------------- + ! Read in ocean forcing data from an existing file + !--------------------------------------------------------------- + write (nu_diag,*) 'ocean mixed layer forcing data file = ', & + trim(sst_file) + write (nu_diag,*) 'FILE YEAR:', myear + call file_year(sst_file,myear) + write (nu_diag,*) 'ocean mixed layer forcing data file = ', & + trim(sst_file) + + endif ! master_task + + if (trim(ocn_data_format) == 'nc') then +#ifdef USE_NETCDF + if (my_task == master_task) then + call ice_open_nc(sst_file, fid) + +! status = nf90_inq_dimid(fid,'nlon',dimid) + status = nf90_inq_dimid(fid,'xt_ocean',dimid) + status = nf90_inquire_dimension(fid,dimid,len=nlon) + +! status = nf90_inq_dimid(fid,'nlat',dimid) + status = nf90_inq_dimid(fid,'yt_ocean',dimid) + status = nf90_inquire_dimension(fid,dimid,len=nlat) + + if( nlon .ne. nx_global ) then + call abort_ice (error_message=subname//'ice: ocn frc file nlon ne nx_global', & + file=__FILE__, line=__LINE__) + endif + if( nlat .ne. ny_global ) then + call abort_ice (error_message=subname//'ice: ocn frc file nlat ne ny_global', & + file=__FILE__, line=__LINE__) + endif + endif ! master_task + ! Read in ocean forcing data for all 12 months + do n=1,ndfld + do m=1,12 + ! Note: netCDF does single to double conversion if necessary + !if (n >= 3 .and. n <= 4) then + ! call ice_read_nc(fid, m, vname(n), work1, debug_forcing, & + ! field_loc_NEcorner, field_type_vector) + if (n == 3 .or. n == 4) then ! 2D currents + nzlev = 1 ! surface currents + !call ice_read_nc_uv(fid, m, nzlev, vname(n), work1, debug_forcing, & + ! field_loc_NEcorner, field_type_vector) + call ice_read_nc(fid, m, vname(n), work1, debug_forcing, & + field_loc_NEcorner, field_type_vector) + else + call ice_read_nc(fid, m, vname(n), work1, debug_forcing, & + field_loc_center, field_type_scalar) + endif + if (n.eq.1) then + ocn_frc_m_access(:,:,:,n,m) = work1(:,:,:) - 273.15 ! Converting from K to C + else + ocn_frc_m_access(:,:,:,n,m) = work1(:,:,:) + endif + enddo ! month loop + enddo ! field loop + + + if (my_task == master_task) call ice_close_nc(fid) +#else + call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined for '//trim(sst_file), & + file=__FILE__, line=__LINE__) +#endif + + + ! do iblk = 1, nblocks + ! do j = 1, ny_block + ! do i = 1, nx_block + ! ocn_frc_m_access(i,j,iblk,1,1) = max(ocn_frc_m_access(i,j,iblk,1,1) ,Tf(i,j,iblk)) ! Making sure sst is not less than freezing temperature + ! ocn_frc_m_access(i,j,iblk,2,1) = max(ocn_frc_m_access(i,j,iblk,2,1), c0) ! Making sure sss is not negative + ! if (ocn_frc_m_access(i,j,iblk,3,1) .gt. c5) then + ! ocn_frc_m_access(i,j,iblk,3,1) = min(ocn_frc_m_access(i,j,iblk,3,1), c5) ! Making sure u-velocity isn't unrealistic + ! elseif (uocn(i,j,iblk) .lt. -c5) then + ! ocn_frc_m_access(i,j,iblk,3,1) = max(ocn_frc_m_access(i,j,iblk,3,1), -c5) ! Making sure v-velocity isn't unrealistic + ! endif + ! if (ocn_frc_m_access(i,j,iblk,4,1) .gt. c5) then + ! ocn_frc_m_access(i,j,iblk,4,1) = c5 ! Making sure v-velocity isn't unrealistic + ! elseif (ocn_frc_m_access(i,j,iblk,4,1) .lt. -c5) then + ! ocn_frc_m_access(i,j,iblk,4,1) = -c5 ! Making sure v-velocity isn't unrealistic + ! endif + ! enddo + ! enddo + ! enddo + +! ND: initialising SSS and SST + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + ! ND: 6/10/23 remove as this causes an SST Jump sst(i,j,iblk) = ocn_frc_m_access(i,j,iblk,1,mmonth) + sss(i,j,iblk) = ocn_frc_m_access(i,j,iblk,2,mmonth) + uocn(i,j,iblk) = ocn_frc_m_access(i,j,iblk,3,mmonth) + vocn(i,j,iblk) = ocn_frc_m_access(i,j,iblk,4,mmonth) + + sss(i,j,iblk) = max(sss(i,j,iblk),c0) + !sst(i,j,iblk) = max(sst(i,j,iblk) - 273.15,Tf(i,j,iblk)) + + ! Making sure current velocities aren't unrealistic + if (uocn(i,j,iblk) .gt. c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .lt. -c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .gt. c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .lt. -c5) uocn(i,j,iblk) = c0 + + + if (runtype.eq.'initial') then ! Freeze/melt is calculated within CICE + ! Initialise to 1.0 just to begin the run + !frzmlt(i,j,iblk) = -1.0_dbl_kind!(-1)*ocn_frc_m_access(i,j,iblk,1,1) + !frzmlt(i,j,iblk) = frzmlt(i,j,iblk) + !if (j .lt. (ny_block/4)) then + ! frzmlt(i,j,iblk) = 1.0_dbl_kind + !else + ! frzmlt(i,j,iblk) = -1.0_dbl_kind + !endif + frzmlt(i,j,iblk) = frzmlt(i,j,iblk) + else + frzmlt(i,j,iblk) = frzmlt(i,j,iblk) + endif + enddo + enddo +enddo + if (debug_forcing) then + if (my_task == master_task) & + write (nu_diag,*) 'ocn_data_access_init' + vmin = global_minval(ocn_frc_m_access(:,:,:,1,1),distrb_info,tmask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,1,1),distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sst',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,2,1),distrb_info,tmask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,2,1),distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sss',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,3,1),distrb_info,umask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,3,1),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'uocn',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,4,1),distrb_info,umask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,4,1),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'vocn',vmin,vmax + vmin = global_minval(frzmlt(:,:,:),distrb_info,umask) + vmax = global_maxval(frzmlt(:,:,:),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'frzmlt',vmin,vmax + + + if (my_task == master_task) & + write (nu_diag,*) 'ocn_data_access_init2' + vmin = global_minval(ocn_frc_m_access(:,:,:,1,2),distrb_info,tmask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,1,2),distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sst',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,2,2),distrb_info,tmask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,2,2),distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sss',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,3,2),distrb_info,umask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,3,2),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'uocn',vmin,vmax + vmin = global_minval(ocn_frc_m_access(:,:,:,4,2),distrb_info,umask) + vmax = global_maxval(ocn_frc_m_access(:,:,:,4,2),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'vocn',vmin,vmax + vmin = global_minval(frzmlt(:,:,:),distrb_info,umask) + vmax = global_maxval(frzmlt(:,:,:),distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'frzmlt',vmin,vmax + + + !if (my_task == master_task) & + ! write (nu_diag,*) 'ocn_data_access_init 12' + !vmin = global_minval(ocn_frc_m_access(:,:,:,1,12),distrb_info,tmask) + !vmax = global_maxval(ocn_frc_m_access(:,:,:,1,12),distrb_info,tmask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst',vmin,vmax + ! vmin = global_minval(ocn_frc_m_access(:,:,:,2,12),distrb_info,tmask) + ! vmax = global_maxval(ocn_frc_m_access(:,:,:,2,12),distrb_info,tmask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sss',vmin,vmax + ! vmin = global_minval(ocn_frc_m_access(:,:,:,3,12),distrb_info,umask) + ! vmax = global_maxval(ocn_frc_m_access(:,:,:,3,12),distrb_info,umask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'uocn',vmin,vmax + ! vmin = global_minval(ocn_frc_m_access(:,:,:,4,12),distrb_info,umask) + ! vmax = global_maxval(ocn_frc_m_access(:,:,:,4,12),distrb_info,umask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'vocn',vmin,vmax + ! vmin = global_minval(frzmlt(:,:,:),distrb_info,umask) + ! vmax = global_maxval(frzmlt(:,:,:),distrb_info,umask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'frzmlt',vmin,vmax + !endif + ! endif + + + !if (my_task == master_task) & + ! write (nu_diag,*) 'ocn_data_access_init point (10,200:,1,1)' + ! vmin = ocn_frc_m_access(10,200,1,1,1) + ! vmax = ocn_frc_m_access(10,200,1,1,1) + !if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst',vmin,vmax + !if (my_task == master_task) & + ! write (nu_diag,*) 'ocn_data_access_init point (10,200:,1,2)' + ! vmin = ocn_frc_m_access(10,200,1,1,2) + ! vmax = ocn_frc_m_access(10,200,1,1,2) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst',vmin,vmax + !if (my_task == master_task) & + ! write (nu_diag,*) 'ocn_data_access_init point (10,200:,1,6)' + ! vmin = ocn_frc_m_access(10,200,1,1,6) + ! vmax = ocn_frc_m_access(10,200,1,1,6) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst',vmin,vmax + endif + + + +! else ! binary format + + !nbits = 64 + !call ice_open (nu_forcing, sst_file, nbits) + + !nrec = 0 + !do n=1,ndfld + ! do m=1,12 + ! nrec = nrec + 1 + ! if (n == 3 .or. n == 4) then ! u-grid + ! call ice_read (nu_forcing, nrec, work1, 'rda8', debug_forcing, & + ! field_loc_NEcorner, field_type_vector) + ! else ! t-grid + ! call ice_read (nu_forcing, nrec, work1, 'rda8', debug_forcing, & + ! field_loc_center, field_type_scalar) + ! endif + ! ocn_frc_m_access(:,:,:,n,m) = work1(:,:,:) + ! enddo ! month loop + !enddo ! field loop + !close (nu_forcing) + + endif + +!echmod - currents cause Fram outflow to be too large +! ocn_frc_m_access(:,:,:,4,:) = c0 +! ocn_frc_m_access(:,:,:,5,:) = c0 +!echmod + +!!!! HYCOM + + !if (my_task == master_task) then + ! write (nu_diag,*)' ' + ! write (nu_diag,*)'Initial ocean forcing file: ',trim(sst_file) + !endif + + !fieldname = 'sss' + !call ice_open_nc (sst_file, fid) + !call ice_read_nc (fid, 1 , fieldname, sss, debug_forcing, & + ! field_loc_center, field_type_scalar) ! t-grid + !call ice_close_nc(fid) + + !call ocn_freezing_temperature + + !sst_file = trim(ocn_data_dir)//'ice.restart.surf.nc' + + !fieldname = 'sst' + !call ice_open_nc (sst_file, fid) + !call ice_read_nc (fid, 1 , fieldname, sst, debug_forcing, & + ! field_loc_center, field_type_scalar) ! t-grid + !call ice_close_nc(fid) + + !fieldname = 'u' + !call ice_open_nc (sst_file, fid) + !call ice_read_nc (fid, 1 , fieldname, uocn, debug_forcing, & + ! field_loc_NEcorner, field_type_vector) ! u-grid + !call ice_close_nc(fid) + + !fieldname = 'v' + !call ice_open_nc (sst_file, fid) + !call ice_read_nc (fid, 1 , fieldname, vocn, debug_forcing, & + ! field_loc_NEcorner, field_type_vector) ! u-grid + !call ice_close_nc(fid) + + ! Make sure sst is not less than freezing temperature Tf + !$OMP PARALLEL DO PRIVATE(iblk,i,j) + !do iblk = 1, nblocks + ! do j = 1, ny_block + ! do i = 1, nx_block + ! sst(i,j,iblk) = max(sst(i,j,iblk) - 273.15,Tf(i,j,iblk)) ! Converting from K to C + ! enddo + ! enddo + !enddo + !$OMP END PARALLEL DO + + + end subroutine ocn_data_access_init + + + +!======================================================================= + + subroutine ocn_data_access(dt) + +! Interpolate monthly ocean data to timestep. +! Restore sst if desired. sst is updated with surface fluxes in ice_ocean.F. + + use ice_blocks, only: nx_block, ny_block + use ice_global_reductions, only: global_minval, global_maxval + use ice_domain, only: nblocks, distrb_info + use ice_domain_size, only: max_blocks + use ice_flux, only: sss, sst, Tf, uocn, vocn, ss_tltx, ss_tlty, & + qdp, hmix + use ice_restart_shared, only: restart + use ice_grid, only: hm, tmask, umask + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + integer (kind=int_kind) :: & + i, j, n, iblk , & + ixm,ixp , & ! record numbers for neighboring months + maxrec , & ! maximum record number + recslot , & ! spline slot for current record + midmonth ! middle day of month + + real (kind=dbl_kind) :: & + vmin, vmax + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + work1 + + character(len=*), parameter :: subname = '(ocn_data_access)' + + if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start' + if (my_task == master_task) & + write (nu_diag,*) 'ocn_data_access' + !------------------------------------------------------------------- + ! monthly data + ! + ! Assume that monthly data values are located in the middle of the + ! month. + !------------------------------------------------------------------- + + midmonth = 15 ! data is given on 15th of every month +! midmonth = fix(p5 * real(daymo(mmonth),kind=dbl_kind)) ! exact middle + + ! Compute record numbers for surrounding months + maxrec = 12 + ixm = mod(mmonth+maxrec-2,maxrec) + 1 + ixp = mod(mmonth, maxrec) + 1 + if (mday >= midmonth) ixm = -99 ! other two points will be used + if (mday < midmonth) ixp = -99 + + if (my_task.eq.master_task) & + write (nu_diag,*) 'ND: ice_forcing ' + if (my_task.eq.master_task) & + write (nu_diag,*) 'mmonth: ',mmonth + if (my_task.eq.master_task) & + write (nu_diag,*) 'ixm, ixp: ',ixm,ixp + + ! Determine whether interpolation will use values 1:2 or 2:3 + ! recslot = 2 means we use values 1:2, with the current value (2) + ! in the second slot + ! recslot = 1 means we use values 2:3, with the current value (2) + ! in the first slot + recslot = 1 ! latter half of month + if (mday < midmonth) recslot = 2 ! first half of month + + ! Find interpolation coefficients + call interp_coeff_monthly (recslot) + + sst_data(:,:,:,:) = c0 + do n = ndfld, 1, -1 + do iblk = 1, nblocks + ! use sst_data arrays as temporary work space until n=1 + if (ixm /= -99) then ! first half of month + do i = 1, nx_block + do j = 1, ny_block + sst_data(i,j,1,iblk) = ocn_frc_m_access(i,j,iblk,n,ixm) + sst_data(i,j,2,iblk) = ocn_frc_m_access(i,j,iblk,n,mmonth) + enddo + enddo + else ! second half of month + sst_data(:,:,1,iblk) = ocn_frc_m_access(:,:,iblk,n,mmonth) + sst_data(:,:,2,iblk) = ocn_frc_m_access(:,:,iblk,n,ixp) + endif + if (debug_forcing) then + if (my_task.eq.master_task) then + !write(nu_diag,*) 'ixm', ixm + !write(nu_diag,*) 'ixp', ixp + !write(nu_diag,*) 'mmonth', mmonth + endif + endif + enddo + + + ! ND: Debugging SST + !if (debug_forcing) then + !if (my_task.eq.master_task) write(nu_diag,*) 'Variable', n + !if (my_task.eq.master_task) write (nu_diag,*) 'Pre interpolation:' + ! vmin = global_minval(sst,distrb_info,tmask) + ! vmax = global_maxval(sst,distrb_info,tmask) + !if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst',vmin,vmax + ! vmin = global_minval(work1,distrb_info,tmask) + ! vmax = global_maxval(work1,distrb_info,tmask) + !if (my_task.eq.master_task) & + ! write (nu_diag,*) 'work1',vmin,vmax + + + ! vmin = global_minval(sst_data(:,:,1,:),distrb_info,tmask) !MAXVAL(sst_data(:,:,1,iblk)) + ! vmax = global_maxval(sst_data(:,:,1,:),distrb_info,tmask) !MINVAL(sst_data(:,:,1,iblk)) + !if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst_data',vmin,vmax + ! vmin = global_minval(sst_data(:,:,2,:),distrb_info,tmask) !MAXVAL(sst_data(:,:,1,iblk)) + ! vmax = global_maxval(sst_data(:,:,2,:),distrb_info,tmask) + !if (my_task.eq.master_task) & + ! write (nu_diag,*) 'sst_data2',vmin,vmax + + + ! vmin = sst_data(10,200,1,1) + ! vmax = sst_data(10,200,1,1) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'ocn_data_access point (10,200,1,1)',vmin,vmax + ! vmin = sst_data(10,200,2,1) + ! vmax = sst_data(10,200,2,1) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'ocn_data_access point (10,200,2,1)',vmin,vmax + !endif !debug_forcing + +! vmin = global_minval(sst_data(:,:,1,iblk),distrb_info,tmask) +! vmax = global_maxval(sst_data(:,:,1,iblk),distrb_info,tmask) +! if (my_task.eq.master_task) & +! write (nu_diag,*) 'sst_data(:,:,1,iblk)',vmin,vmax +! vmin = global_minval(sst_data(:,:,2,iblk),distrb_info,tmask) +! vmax = global_maxval(sst_data(:,:,2,iblk),distrb_info,tmask) +! if (my_task.eq.master_task) & +! write (nu_diag,*) 'sst_data(:,:,2,iblk)',vmin,vmax + + + call interpolate_data (sst_data,work1) + + ! if (debug_forcing) then + !vmin = work1(10,200,iblk) + !vmax = work1(10,200,iblk) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'work1 interp: work1 ',vmin,vmax + + ! vmin = global_minval(work1(:,:,:),distrb_info,tmask) + ! vmax = global_maxval(work1(:,:,:),distrb_info,tmask) + ! if (my_task.eq.master_task) & + ! write (nu_diag,*) 'POST interp: work1 ',vmin,vmax + ! endif + + ! masking by hm is necessary due to NaNs in the data file + do j = 1, ny_block + do i = 1, nx_block + if (n == 2) sss (i,j,:) = c0 + if (n == 3) uocn (i,j,:) = c0 + if (n == 4) vocn (i,j,:) = c0 + do iblk = 1, nblocks + if (hm(i,j,iblk) == c1) then + if (n == 2) sss (i,j,iblk) = work1(i,j,iblk) + if (n == 3) uocn (i,j,iblk) = work1(i,j,iblk) + if (n == 4) vocn (i,j,iblk) = work1(i,j,iblk) + endif + enddo + enddo + enddo + enddo + + + + do j = 1, ny_block + do i = 1, nx_block + sss (i,j,:) = max (sss(i,j,:), c0) + enddo + enddo + + call ocn_freezing_temperature + + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + ! work1(i,j,iblk) = max(work1(i,j,iblk) ,-c5) ! Making sure sst is not less than freezing temperature + work1(i,j,iblk) = max(sst(i,j,iblk), Tf(i,j,iblk)) + enddo + enddo + enddo + + + + if (restore_ocn) then + do j = 1, ny_block + do i = 1, nx_block + sst(i,j,:) = sst(i,j,:) + (work1(i,j,:)-sst(i,j,:))*dt/trest + enddo + enddo +! else sst is only updated in ice_ocean.F + endif + + ! initialize sst properly on first step + if (istep1 <= 1 .and. .not. (restart)) then + call interpolate_data (sst_data,sst) + !$OMP PARALLEL DO PRIVATE(iblk,i,j) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + if (hm(i,j,iblk) == c1) then + sst(i,j,iblk) = max (sst(i,j,iblk), Tf(i,j,iblk)) + else ! ND: It is land! + sst(i,j,iblk) = c0 + sss(i,j,iblk) = c0 + uocn(i,j,iblk) = c0 + vocn(i,j,iblk) = c0 + endif + !if (umask(i,j,iblk) .eqv. .false.) then + ! ! ND: Setting ocean velocities on land to 0 + ! uocn(i,j,iblk) = c0 + ! vocn(i,j,iblk) = c0 + !endif + ! Turning off velocities in the NH + !if (j.gt.(ny_block/2)) uocn(i,j,iblk) = c0 + !if (j.gt.(ny_block/2)) vocn(i,j,iblk) = c0 + ! Making sure current velocities aren't unrealistic + if (uocn(i,j,iblk) .gt. c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .lt. -c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .gt. c5) uocn(i,j,iblk) = c0 + if (uocn(i,j,iblk) .lt. -c5) uocn(i,j,iblk) = c0 + enddo + enddo + enddo + !else + ! do iblk = 1, nblocks + ! do j = 1, ny_block + ! do i = 1, nx_block + ! if (hm(i,j,iblk) == c1) then + ! sst(i,j,iblk) = max (sst(i,j,iblk), Tf(i,j,iblk)) + ! else ! ND: It is land! + ! sst(i,j,iblk) = c0 + ! sss(i,j,iblk) = c0 + ! uocn(i,j,iblk) = c0 + ! vocn(i,j,iblk) = c0 + ! endif + ! if (umask(i,j,iblk) .eqv. .false.) then + ! ! ND: Setting ocean velocities on land to 0 + ! uocn(i,j,iblk) = c0 + ! vocn(i,j,iblk) = c0 + ! endif + ! if (j.gt.(ny_block/2)) uocn(i,j,iblk) = c0 + ! if (j.gt.(ny_block/2)) vocn(i,j,iblk) = c0 + ! enddo + ! enddo + !enddo + !$OMP END PARALLEL DO + endif + + if (debug_forcing) then + if (my_task == master_task) & + write (nu_diag,*) 'ocn_data_access' + vmin = global_minval(sst,distrb_info,tmask) + vmax = global_maxval(sst,distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sst',vmin,vmax + vmin = global_minval(sss,distrb_info,tmask) + vmax = global_maxval(sss,distrb_info,tmask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'sss',vmin,vmax + vmin = global_minval(uocn,distrb_info,umask) + vmax = global_maxval(uocn,distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'uocn',vmin,vmax + vmin = global_minval(vocn,distrb_info,umask) + vmax = global_maxval(vocn,distrb_info,umask) + if (my_task.eq.master_task) & + write (nu_diag,*) 'vocn',vmin,vmax + endif + + end subroutine ocn_data_access + +!======================================================================= + + + + +end module ice_forcing + +!======================================================================= +>>>>>>> Stashed changes diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90 index d65cf52d3..54d83d663 100644 --- a/cicecore/cicedynB/general/ice_step_mod.F90 +++ b/cicecore/cicedynB/general/ice_step_mod.F90 @@ -543,7 +543,7 @@ subroutine step_therm2 (dt, iblk) use ice_flux_bgc, only: flux_bio, faero_ocn, & fiso_ocn, HDO_ocn, H2_16O_ocn, H2_18O_ocn use ice_grid, only: tmask - use ice_state, only: aice, aicen, aice0, trcr_depend, & + use ice_state, only: aice, vice, aicen, aice0, trcr_depend, & aicen_init, vicen_init, trcrn, vicen, vsnon, & trcr_base, n_trcr_strata, nt_strata @@ -572,7 +572,74 @@ subroutine step_therm2 (dt, iblk) character(len=*), parameter :: subname = '(step_therm2)' +<<<<<<< Updated upstream call icepack_query_parameters(z_tracers_out=z_tracers,solve_zsal_out=solve_zsal) +======= + ! Noah Day WIM------------------------------------------------------------ + real(kind=dbl_kind), dimension(nfreq) :: alpha_coeff, wave_spectrum_in, wave_spectrum_out ! attenuation coefficient + real(kind=dbl_kind) :: hice_init, Length_cell, aice_in!conc_coeff, Dav, Lcell, m0, m2, ice_thick + real (kind=dbl_kind), dimension(nx_block,ny_block) :: worka ! work for calculating mean FSD + real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,nblocks) :: wave_spec_blk ! wave spectrum in m^2s/rad + + real(kind=dbl_kind), dimension(nfreq) :: om ! freqs, rad/s + real(kind=dbl_kind) :: fmin, fmax + ! freq min/max, 1/s + real(kind=dbl_kind) :: om1, om2 + ! ang freqs, rad/s + real(kind=dbl_kind) :: om_0 + + real(kind=dbl_kind) :: loc_swh ! dummy swh holder + + real (kind=dbl_kind), dimension(1) :: & + temp_loc + + + integer (kind=int_kind) :: ii, jj, kk, n, k ! index + integer :: lp_i, lp_j ! counter + integer (kind=int_kind) :: nthh, tmtt, idll + + ! Wave-ice variables: + real (kind=dbl_kind), dimension(nx_block,ny_block) :: & + ov_vol, ov_conc + + integer (kind=int_kind), save :: & + icells ! number of cells with aicen > puny + ! wavcells ! LB (Oct `13): number of cells to apply wave-ice code in + + ! integer, parameter :: & + ! nedge=5 ! LB (Oct `13): ice-edge + + integer (kind=4) :: & + wavemask_dyn ! LB (Sept `14): dynamic ice edge/wavemask + + integer (kind=int_kind), dimension(nx_block) :: & + wavemask_dyn_vec, mn_lat_vec, dumlatloc_vec !Noah Day, dynamic ice edge/wavemask vector (for each block) + + integer (kind=int_kind), dimension(nx_block*ny_block) :: & + indxi, indxj ! indirect indices for cells with aicen > puny + + ! LB: added 20.05.15 for WW3 data integration + integer, dimension(1) :: dumlatloc + real (kind=dbl_kind) :: mn_lat + + integer (kind=int_kind) :: floe_count ! whether we have memory or not (0=none, 1=yes) + + !real (kind=dbl_kind) :: & + !uvel_center, & ! cell-centered velocity, x component (m/s) + !vvel_center, & ! cell-centered velocity, y component (m/s) + !puny ! a very small number + !----------------------------------------------------------------- + ! ice state at start of time step, saved for later in the step + !----------------------------------------------------------------- + + !real (kind=dbl_kind), dimension(nx_block,ny_block,nblocks) :: & + ! aice_init, & ! initial concentration of ice, for diagnostics + ! floediam_init ! initial floe diameter, for diagnostics (LB 08.07.14) + + ! ------------------------------------------------------------------------ + + call icepack_query_parameters(z_tracers_out=z_tracers,solve_zsal_out=solve_zsal, puny_out=puny) +>>>>>>> Stashed changes call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) @@ -591,15 +658,396 @@ subroutine step_therm2 (dt, iblk) ihi = this_block%ihi jlo = this_block%jlo jhi = this_block%jhi +<<<<<<< Updated upstream +======= + if (cmt.ne.0) then + write(nu_diag,*) ' ---- > Calculating for block: ', iblk + write(nu_diag,*) ' ilo is: ', ilo + write(nu_diag,*) ' ihi is: ', ihi + write(nu_diag,*) ' jlo is: ', jlo + write(nu_diag,*) ' jhi is: ', jhi + write(nu_diag,*) ' ----------------------------------------> ' + write(nu_diag,*) ' WIM .eq. ', WIM + end if + +!------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! TURNING ON WAVES IN ICE MODULE +!------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- + call ice_timer_start(timer_wim) + + if (WIM.eq.1) then + !----------------------------------------------------------------- + ! Save the ice area passed to the coupler (so that history fields + ! can be made consistent with coupler fields). + ! Save the initial ice area and volume in each category. + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! MIZ FLOE SIZE DISTRIBUTION: BASED ON OCEAN WAVES - Noah Day WIM + !----------------------------------------------------------------- + + ! find cells with aice <= puny + !icells = 0 + wavemask_dyn = 0 !216 ! ND: 06-03-23 was 0, changing to 57.5 S + do jj = jlo, jhi + do ii = ilo, ihi + if (aice(ii,jj,iblk).gt.puny) then + !icells = icells + 1 + !indxi(icells) = ii + !indxj(icells) = jj + if ((jj.gt.wavemask_dyn).and.(jj.lt.ny_block/2)) then + wavemask_dyn = jj + endif + endif + enddo ! ii + enddo ! jj + + if (cmt.ne.0) then + write(nu_diag,*) 'indxi: ', ii + write(nu_diag,*) ' indxj: ', jj + write(nu_diag,*) 'wave mask is: ', wavemask_dyn + end if !cmt + + if (max_floediam.eq.300.0_dbl_kind) then !!! Wave-ice interaction code ON - Noah Day commenting out 26/03/22 + + if (cmt.ne.0) then + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) 'LB: CICE_RunMod -> Running wave-ice interaction code:', & + ' istep=', istep !, 'dt=', dt + endif ! cmt + + ! Calculate mean latitude + mn_lat = sum(TLAT(:,wavemask_dyn,1))/size(TLAT(:,wavemask_dyn,1)) + if (cmt.ne.0) write(nu_diag,*) ' mean latitude: ', mn_lat, mn_lat*c180/pi + + if (WAVE_METH.eq.1) then + ! Finding the lateral degree where to propagate from + dumlatloc=minloc(abs(ww3_lat-c180*mn_lat/pi),dim=2) ! radians to degrees for mean lat + if (cmt.ne.0) write(nu_diag,*)' dumlatloc: ', dumlatloc + + if (WIM_LONG.ne.0) then + !write(nu_diag,*) ' Longitudinal WIM Starting.....' + ! Calculate the wavemask for each longitude + wavemask_dyn_vec(:) = 216 ! ND: 06-03-23 was 0, changing to 57.5 S + do jj = jlo, jhi + do ii = ilo, ihi + if (aice(ii,jj,iblk).gt.puny) then + !icells = icells + 1 + !indxi(icells) = ii + !indxj(icells) = jj + if ((jj.lt.ny_block/2).and.(jj.gt.0)) then + wavemask_dyn_vec(ii) = jj + endif + endif + end do + end do + +! do ii = ilo, ihi +! do jj = jlo, jhi! Start the latitude at the bottom of the grid +! if (tmask(ii,jj,iblk)) then +! icells = icells + 1 +! indxi(icells) = ii +! indxj(icells) = jj +! if (aice(ii,jj,iblk).lt.puny) then +! if ((jj.lt.ny_block/2).and.(jj.gt.0)) wavemask_dyn_vec(ii) = jj +! exit +! endif +! endif +! end do +! end do +! write(nu_diag,*) 'iblk', iblk +! write(nu_diag,*) 'wavemask_dyn_vec', wavemask_dyn_vec + + ! ND: FIX THIS FOR DIFFERENT NUMBER OF HALO CELLS + wavemask_dyn_vec(1) = wavemask_dyn_vec(2) + wavemask_dyn_vec(nx_block) = wavemask_dyn_vec(nx_block-1) + ! Initialise the dummy latitude location vector + dumlatloc_vec(:) = 0 + + do ii = ilo, ihi + temp_loc = minloc(abs(ww3_lat-c180*mn_lat_vec(ii)/pi),dim=2) + dumlatloc_vec(ii) = temp_loc(1) + !write(nu_diag,*)' dumlatloc: ', minloc(abs(ww3_lat-c180*mn_lat_vec(ii)/pi),dim=2) + end do + + ! Initialise wave spectrum per longitude (equivalent to cell for 1 degree grid) + call init_wave_spec_long(wavemask_dyn_vec,ww3_swh(:,:), & + ww3_fp(:,:),ww3_dir(:,:),size(ww3_lon),size(ww3_lat),iblk) + + if (OVERWRITE_DIRS.eq.1) then ! Force waves to the south + !print*, 'OVERWRITING directions: wavemask_dyn=', wavemask_dyn + do i=1,nx_block + mwd(i,wavemask_dyn_vec(i),iblk) = pi + enddo + endif ! IF OVERWRITE_DIRS + + else ! WIM_LONG: Run blockwise + !write(nu_diag,*) ' Blockwise WIM.....' + ! Initialise wave spectrum + call init_wave_spec(wavemask_dyn,ww3_swh(:,dumlatloc(1)), & + ww3_fp(:,dumlatloc(1)),ww3_dir(:,dumlatloc(1)),size(ww3_lon),iblk) + if (OVERWRITE_DIRS.eq.1) then ! Force waves to the south + !print*, 'OVERWRITING directions: wavemask_dyn=', wavemask_dyn + do i=1,nx_block + mwd(i,wavemask_dyn,iblk) = pi + enddo + endif ! IF OVERWRITE_DIRS + endif ! WIM_LONG + !print*, 'passed overwriting directions: mwd(1,wavemask_dyn,iblk)=', mwd(1,wavemask_dyn,iblk) + else ! WAVE_METH + call init_wave_spec_usr(wavemask_dyn) + ! write(nu_diag,*) ' -> Calling init_wave_spec_usr' + endif ! WAVE_METH + + if (1.eq.0) then ! reinitialise floe diameters (no memory) + !if (idate.eq.idate0) then!.and.timesecs.eq.sec_init) then + !if (istep.eq.1) then + if (cmt.ne.0) write(nu_diag,*) ' -> Reinitialising floe sizes' + call init_floe_0 + else ! initialise ifd using ifloe tracer values + if (cmt.ne.0) write(nu_diag,*) ' -> Remembering floe sizes' + ! Noah Day ifd(:,:,iblk) = trcrn(:,:,nt_fsd,1,iblk) + ! Noah Day 20/3/22 START --------------------------------------- + do i = 1,nx_block + do j = 1,ny_block + if (aice(i,j,iblk).gt.puny) then + !write(nu_diag,*) 'For cell i,j:',i,j + !write(nu_diag,*) 'trcrn(:,:,nt_fsd,1,iblk)', trcrn(i,j,nt_fsd,1,iblk) + worka(i,j) = c0 + ! Calculate the representative radius (Roach et al. 2018) + do k = 1, nfsd + do n = 1, ncat + worka(i,j) = worka(i,j) & + + (trcrn(i,j,nt_fsd+k-1,n,iblk) * floe_rad_c(k) & + * aicen(i,j,n,iblk)/aice(i,j,iblk)) + end do + end do + + ! Feed in the diameter to WIM (Bennetts et al. 2017) + ifd(i,j,iblk) = c2*worka(i,j) + else ! Otherwise there is no ice + ifd(i,j,iblk) = c0 + end if + end do + end do + ! Noah Day 20/3/22 END --------------------------------------- + endif ! date + + ! Initialise overall volume & overall concentration (wrt thickness categories) + do j = 1, ny_block + do i = 1, nx_block + ov_conc(i,j) = c0 + ov_vol(i,j) = c0 + enddo + enddo + + do j = wavemask,2,-1 + do i = 1, nx_block + if (tmask(i,j,iblk)) then + do n = 1, ncat + ov_conc(i,j) = ov_conc(i,j) + aicen(i,j,n,iblk) + ov_vol(i,j) = ov_vol(i,j) + vicen(i,j,n,iblk) + enddo + endif + enddo + enddo + + + do j = 1, ny_block + do i = 1, nx_block + ov_conc(i,j) = aice(i,j,iblk) + ov_vol(i,j) = vice(i,j,iblk) + enddo + enddo + + ! Convert from m^2s to m^2s/rad + do j = 1, ny_block + do i = 1, nx_block + do k = 1, nfreq + wave_spec_blk(i,j,k,iblk) = wave_spectrum(i,j,k,iblk)/(2*pi) + enddo + enddo + enddo + + if (WIM_LONG.eq.1) then ! Propagate waves per cell + !write(nu_diag,*) 'Calling increment_floe_long..............' + call increment_floe_long (nx_block, ny_block, & ! nx_block, ny_block + dt, & ! dt + tmask(:,:,iblk), & ! tmask + HTE(:,:,iblk), & ! Lcell + swh(:,:,iblk), & ! loc_swh + ppd(:,:,iblk), & ! loc_ppd + mwd(:,:,iblk), & ! loc_mwd + ifd(:,:,iblk), & ! ifloe + ov_conc, ov_vol, & ! afice, vfice + wavemask_dyn, & ! dum_wavemask + wavemask_dyn_vec, & ! dum_wavemask + wave_spec_blk(:,:,:,iblk)) ! wave spectrum + else if (WIM_LONG.eq.0) then ! Propagate waves on a block wise basis + call increment_floe (nx_block, ny_block, & ! nx_block, ny_block + dt, & ! dt + tmask(:,:,iblk), & ! tmask + HTE(:,:,iblk), & ! Lcell + swh(:,:,iblk), & ! loc_swh + ppd(:,:,iblk), & ! loc_ppd + mwd(:,:,iblk), & ! loc_mwd + ifd(:,:,iblk), & ! ifloe + ov_conc, ov_vol, & ! afice, vfice + wavemask_dyn, & ! dum_wavemask + wave_spec_blk(:,:,:,iblk)) ! wave spectrum + !write(nu_diag,*) ' Called increment floe' + endif ! WIM_LONG + ! write(nu_diag,*) '----------------------------------------------------' + ! write(nu_diag,*) '------------WAVE PROP DONE----------------------' + ! Convert from m^2s/rad to m^2s + do j = 1, ny_block + do i = 1, nx_block + do k = 1, nfreq + wave_spectrum(i,j,k,iblk) = wave_spec_blk(i,j,k,iblk)*(2*pi) + ! if (k.gt.1.0) write(nu_diag,*) ' Cell with waves has coords:', i,j + enddo + enddo + enddo + + else !!! Wave-ice interaction code OFF + + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) '----------------------------------------------------' + write(nu_diag,*) 'LB: CICE_RunMod > Control test -> setting floe diameters = ', & + max_floediam !c300 + endif ! ENDIF c300 + + ! Noah Day WIM 25/10/21 + ! Creating a wave spectrum in the icepack using the attenuated significant wave height and peak period + fmin = 1d0/1000d0 ! minimum frequency, 0.001 Hz + fmax = 1d0/1d0 ! maximum frequency, 1 Hz + om1=2*pi*fmin ! angular frequency, rad/s + om2=2*pi*fmax ! angular frequency, rad/s + om_0 = (om2 - om1)/(nfreq-1) ! nfreq must match nw in WIM + ! Importing om values + + do lp_i=1,nfreq + om(lp_i) = om1 + (lp_i-1)*om_0 + wavefreq(lp_i) = om(lp_i)/(c2*pi) ! Converting from angular frequency to Hz + end do + + dwavefreq(1) = wavefreq(1) + do lp_i=2,nfreq + dwavefreq(lp_i) = wavefreq(lp_i) - wavefreq(lp_i-1) + end do + + + ! TESTING + ! do j = 1, ny_block + ! do i = 1, nx_block + ! write(nu_diag,*) 'Freq SWH:', c4* SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:))) + ! write(nu_diag,*) 'Freq SWH 2:', c4* SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:)*c2*pi)) + ! write(nu_diag,*) 'Ang Freq SWH:', c4* SQRT(SUM(wave_spec_blk(i,j,:,iblk)*om_0)) + ! write(nu_diag,*) 'SWH:', swh(i,j,iblk) + ! if (k.gt.1.0) write(nu_diag,*) ' Cell with waves has coords:', i,j + ! enddo + ! enddo + + + !dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind)) + + !------------------------------------------------------------------------------- + ! TURNING OFF WAVES IN ICE MODULE + call ice_timer_stop(timer_wim) + endif ! WIM + !------------------------------------------------------------------------------- +>>>>>>> Stashed changes do j = jlo, jhi do i = ilo, ihi if (tmask(i,j,iblk)) then +<<<<<<< Updated upstream ! significant wave height for FSD if (tr_fsd) & wave_sig_ht(i,j,iblk) = c4*SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:))) +======= + ! significant wave height for FSD + if (tr_fsd) then + if (WIM.eq.1) then ! Noah Day WIM, if WIM is off: calculate SWH from dummy data + if (WAVE_METH.eq.0) then ! Use dummy data from CICE code + ! set for 25 frequencies + !nfreq = 25 ! ND: previously set to 25 + wave_spectrum(i,j,:,iblk) = c0 + + ! FOR TESTING ONLY - do not use for actual runs!! + wave_spectrum(i,j,1,iblk) = 0.00015429197810590267 + wave_spectrum(i,j,2,iblk) = 0.002913531381636858 + wave_spectrum(i,j,3,iblk) = 0.02312942035496235 + wave_spectrum(i,j,4,iblk) = 0.07201970368623734 + wave_spectrum(i,j,5,iblk) = 0.06766948103904724 + wave_spectrum(i,j,6,iblk) = 0.005527883302420378 + wave_spectrum(i,j,7,iblk) = 3.326293881400488e-05 + wave_spectrum(i,j,8,iblk) = 6.815936703929992e-10 + wave_spectrum(i,j,9,iblk) = 2.419401186610744e-20 + ! hardwired for wave coupling with NIWA version of Wavewatch + ! From Wavewatch, f(n+1) = C*f(n) where C is a constant set by the user + ! These freq are for C = 1.1 + wavefreq = (/0.04118, 0.045298, 0.0498278, 0.05481058, 0.06029164, & + 0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029, & + 0.10681032, 0.11749136, 0.1292405, 0.14216454, 0.15638101, & + 0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317, & + 0.27703848, 0.30474234, 0.33521661, 0.36873826, 0.40561208/) + + ! boundaries of bin n are at f(n)*sqrt(1/C) and f(n)*sqrt(C) + dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind)) + + wave_sig_ht(i,j,iblk) = c4*SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:))) + else ! WAVE_METH: Data from WIM wave propagation + ! Update SWH and PPD + wave_sig_ht(i,j,iblk) = swh(i,j,iblk) + peak_period(i,j,iblk) = ppd(i,j,iblk) + mean_wave_dir(i,j,iblk) = mwd(i,j,iblk) + endif ! WAVE_METH + else ! WIM == 0: No waves + !write(nu_diag,*) 'ND: NO WAVES BEING FORCED' + wave_sig_ht(i,j,iblk) = c0 + peak_period(i,j,iblk) = c0 + mean_wave_dir(i,j,iblk) = c0 + endif ! WIM + endif ! tr_fsd + + ! Using Bretschneider to create a wave spectrum + if (wave_sig_ht(i,j,iblk).gt.puny.and.peak_period(i,j,iblk).gt.puny) then + if (cmt.ne.0) then + write(nu_diag,*) 'FOR CELL (i,j,iblk): ', i,j,iblk + write(nu_diag,*) ' ---- > SWH(i,j,iblk) before: ', wave_sig_ht(i,j,iblk) + write(nu_diag,*) ' ---- > wave_spectrum SHAPE: ', SHAPE(wave_spectrum) + write(nu_diag,*) ' ---- > dwavefreq: ', dwavefreq + write(nu_diag,*) ' ---- > wave_spectrum(i,j,:,iblk): ', wave_spectrum(i,j,:,iblk) + end if ! cmt + !!do lp_i=1,nfreq ! calculate the wave spectrum for each cell given peak period and significant wave height + !! wave_spectrum(i,j,lp_i,iblk) = SDF_Bretschneider(om(lp_i),0,swh(i,j,iblk),peak_period(i,j,iblk)) + !!end do + !loc_swh = fn_SpecMoment(swh(i,j,:,iblk),nfreq,1,om,mean_wave_dir(i,j,iblk),0,nu_diag) + ! c4*(loc_swh**(p5)) ! recalculating the SWH from the Spectrum at that point + if (cmt.ne.0) then + !write(nu_diag,*) ' ---- > SWH(i,j,iblk) after : ', wave_sig_ht(i,j,iblk) + !write(nu_diag,*) ' ---- > wave_spectrum(i,j,:,iblk): ', wave_spectrum(i,j,:,iblk) + end if ! cmt + else ! SWH approx 0 + wave_spectrum(i,j,:,iblk) = c0 + !peak_period(i,j,:,iblk) = c0 + !mean_wave_dir(i,j,iblk) = pi + end if ! SWH + Tp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!! END WAVE-ICE CODE !!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! END NOAH DAY ----------------------------------------------------------------- +! Back to CICE 6 code + !write(nu_diag,*) 'WIM COMPLETE !!!!!!!!!!!!!!!!! ' +>>>>>>> Stashed changes call icepack_step_therm2(dt=dt, ncat=ncat, & nltrcr=nltrcr, nilyr=nilyr, nslyr=nslyr, nblyr=nblyr, & diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 9ce22a6f4..291e80d8e 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -262,10 +262,17 @@ subroutine init_grid1 trim(grid_type) == 'regional' ) then if (trim(grid_format) == 'nc') then +<<<<<<< Updated upstream call ice_open_nc(grid_file,fid_grid) call ice_open_nc(kmt_file,fid_kmt) +======= + !write(nu_diag,*) 'ND: ice_grid.F90 pre grid and kmt' + call ice_open_nc(grid_file,fid_grid) + call ice_open_nc(kmt_file,fid_kmt) + !write(nu_diag,*) 'ND: ice_grid.F90 opened grid and kmt' +>>>>>>> Stashed changes fieldname='ulat' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.) fieldname='kmt' diff --git a/cicecore/cicedynB/infrastructure/ice_read_write.F90 b/cicecore/cicedynB/infrastructure/ice_read_write.F90 index d902c62f8..81acc4019 100644 --- a/cicecore/cicedynB/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedynB/infrastructure/ice_read_write.F90 @@ -1164,7 +1164,9 @@ subroutine ice_read_nc_xy(fid, nrec, varname, work, diag, & !------------------------------------------------------------- status = nf90_inq_varid(fid, trim(varname), varid) - + !write(nu_diag,*) 'ND: varname ', trim(varname) + !write(nu_diag,*) 'ND: status ', status + !write(nu_diag,*) 'ND: fid ', fid if (status /= nf90_noerr) then call abort_ice (subname//'ERROR: Cannot find variable '//trim(varname) ) endif @@ -2137,11 +2139,29 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) count=(/nx_global,ny_global,1/) ) endif endif ! my_task = master_task +<<<<<<< Updated upstream +======= +#ifndef ORCA_GRID + status = nf90_get_var( fid, varid, work_g, & + start=(/1,1,nrec/), & + count=(/nx_global,ny_global,1/) ) +#else + status = nf90_get_var( fid, varid, work_g3)!, & + work_g = work_g3!(2:nx_global+1,1:ny_global) +#endif + !endif ! my_task = master_task + + +>>>>>>> Stashed changes !------------------------------------------------------------------- ! optional diagnostics !------------------------------------------------------------------- +<<<<<<< Updated upstream +======= + !write(nu_diag,*) 'ND: optional diagnostics ' +>>>>>>> Stashed changes if (my_task == master_task .and. diag) then ! write(nu_diag,*) & ! 'ice_read_global_nc, fid= ',fid, ', nrec = ',nrec, & diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 index 1a5681b38..da87bf0b0 100644 --- a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 +++ b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 @@ -19,7 +19,7 @@ module ice_restart_driver use ice_kinds_mod use ice_arrays_column, only: oceanmixed_ice - use ice_constants, only: c0, c1, p5, & + use ice_constants, only: c0, c1, p5, p001, & ! ND: adding p001 field_loc_center, field_loc_NEcorner, & field_type_scalar, field_type_vector use ice_restart_shared, only: restart_dir, pointer_file, & @@ -206,7 +206,7 @@ subroutine restartfile (ice_ic) strocnxT, strocnyT, sst, frzmlt, iceumask, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & - stress12_1, stress12_2, stress12_3, stress12_4 + stress12_1, stress12_2, stress12_3, stress12_4, Tf ! ND: adding tf use ice_flux, only: coszen use ice_grid, only: tmask, grid_type use ice_state, only: trcr_depend, aice, vice, vsno, trcr, & @@ -421,11 +421,19 @@ subroutine restartfile (ice_ic) if (my_task == master_task) & write(nu_diag,*) 'min/max sst, frzmlt' - call read_restart_field(nu_restart,0,sst,'ruf8', & 'sst',1,diag,field_loc_center, field_type_scalar) call read_restart_field(nu_restart,0,frzmlt,'ruf8', & 'frzmlt',1,diag,field_loc_center, field_type_scalar) + ! ND: constant initialisation when using ACCESS-OM2 restart files + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + !frzmlt(i,j,iblk) = -p001 ! ND: Initialise some freeze melt potential to kickstart the ocean (breaks with 0) + sst(i,j,iblk) = max (sst(i,j,iblk), Tf(i,j,iblk)) ! ND: Initialise freezing everywhere + enddo + enddo + enddo endif !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 index 9c6b30ee1..8e5b5b5c0 100644 --- a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 @@ -44,7 +44,7 @@ module ice_history_write subroutine ice_write_hist (ns) use ice_kinds_mod - use ice_arrays_column, only: hin_max, floe_rad_c + use ice_arrays_column, only: hin_max, floe_rad_c, floe_binwidth use ice_blocks, only: nx_block, ny_block use ice_broadcast, only: broadcast_scalar use ice_calendar, only: msec, timesecs, idate, idate0, write_ic, & @@ -75,7 +75,7 @@ subroutine ice_write_hist (ns) integer (kind=int_kind) :: i,k,ic,n,nn, & ncid,status,imtid,jmtid,kmtidi,kmtids,kmtidb, cmtid,timid,varid, & - nvertexid,ivertex,kmtida,iflag, fmtid + nvertexid,ivertex,kmtida,iflag, fmtid, wdhid ! ND: adding wdhid integer (kind=int_kind), dimension(3) :: dimid integer (kind=int_kind), dimension(4) :: dimidz integer (kind=int_kind), dimension(5) :: dimidcz @@ -293,6 +293,9 @@ subroutine ice_write_hist (ns) var_nz(4) = coord_attributes('VGRDb', 'vertical ice-bio levels', '1') var_nz(5) = coord_attributes('VGRDa', 'vertical snow-ice-bio levels', '1') var_nz(6) = coord_attributes('NFSD', 'category floe size (center)', 'm') + !var_nz(7) = coord_attributes('WFSD', 'category floe size (width)', 'm') + + !----------------------------------------------------------------- ! define information for optional time-invariant variables @@ -388,13 +391,14 @@ subroutine ice_write_hist (ns) endif enddo - ! Extra dimensions (NCAT, NZILYR, NZSLYR, NZBLYR, NZALYR, NFSD) + ! Extra dimensions (NCAT, NZILYR, NZSLYR, NZBLYR, NZALYR, NFSD, WFSD) dimidex(1)=cmtid dimidex(2)=kmtidi dimidex(3)=kmtids dimidex(4)=kmtidb dimidex(5)=kmtida dimidex(6)=fmtid + !dimidex(7)=wdhid ! ND: adding wdhid for floe bin width do i = 1, nvarz if (igrdz(i)) then @@ -1170,7 +1174,9 @@ subroutine ice_write_hist (ns) CASE ('NCAT') status = nf90_put_var(ncid,varid,hin_max(1:ncat_hist)) CASE ('NFSD') - status = nf90_put_var(ncid,varid,floe_rad_c(1:nfsd_hist)) + status = nf90_put_var(ncid,varid,floe_rad_c(1:nfsd_hist)) + !CASE ('WFSD') ! ND: floe binwidth + ! status = nf90_put_var(ncid,varid,floe_binwidth(1:nfsd_hist)) CASE ('VGRDi') ! index - needed for Met Office analysis code status = nf90_put_var(ncid,varid,(/(k, k=1,nzilyr)/)) CASE ('VGRDs') ! index - needed for Met Office analysis code diff --git a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 index 08059435f..8a77012b6 100644 --- a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 +++ b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 @@ -24,6 +24,23 @@ module CICE_RunMod use icepack_intfc, only: icepack_max_iso, icepack_max_aero use icepack_intfc, only: icepack_query_parameters use icepack_intfc, only: icepack_query_tracer_flags, icepack_query_tracer_sizes +<<<<<<< Updated upstream +======= +! Noah Day WIM, ---------------------------------------------------------------- + use m_prams_waveice, only: waveicedatadir, fname_ww3, WAVE_METH, ww3_lat, ww3_lon, & + ww3_dir, ww3_tm, ww3_swh, ww3_fp, ATTEN_METH, ATTEN_MODEL, attn_fac, do_coupled, & + OVERWRITE_DIRS, ww3_dir_full, ww3_swh_full, ww3_fp_full, nww3_dt, WIM + use netcdf + use ice_constants, only: eps1, eps3 + use ice_fileunits + use ice_read_write + use ice_restart_shared, only: lenstr, restart_dir, restart_file, & + pointer_file, runtype + use ice_communicate, only: my_task, master_task + use ice_exit, only: abort_ice + use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global ! Noah Day WIM +! ----------------------------------------------------------------------------- +>>>>>>> Stashed changes implicit none private @@ -45,7 +62,11 @@ subroutine CICE_Run use ice_calendar, only: istep, istep1, dt, stop_now, advance_timestep use ice_forcing, only: get_forcing_atmo, get_forcing_ocn, & +<<<<<<< Updated upstream get_wave_spec +======= + get_wave_spec, init_wave_block, check ! Noah Day WIM adding init_wave_block and scheck +>>>>>>> Stashed changes use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & fiso_default, faero_default use ice_flux, only: init_flux_atm, init_flux_ocn @@ -77,11 +98,94 @@ subroutine CICE_Run ! timestep loop !-------------------------------------------------------------------- +<<<<<<< Updated upstream timeLoop: do -#endif +======= +! Noah Day WIM, ---------------------------------------------------------------- +if (WIM.eq.1) then + ! Initialise the indexes for WW3 reading. + nww3 = 24*(mday-1)! Noah Day 9/2 temporary fix this 0 + nmth = mmonth + nyr = myear + write(nu_diag,*) ' Starting date is: ' + write(nu_diag,*) 'mmonth: ', mmonth + write(nu_diag,*) 'mday: ', mday + write(nu_diag,*) 'myear: ', myear + write(nu_diag,*) 'idate: ', idate + if (WAVE_METH.eq.1) then + call sub_WW3_dataread(nmth,N_tm,N_lat,N_lon,nyr,mday) + allocate(ww3_swh(N_lon,N_lat)) + allocate(ww3_fp(N_lon,N_lat)) + allocate(ww3_dir(N_lon,N_lat)) + endif + ! Distribute the wave data across the blocks + !call init_wave_block(ww3_swh,ww3_fp,ww3_dir,N_lon,N_lat) + !nww3=1-nww3_dt + !nww3=1-nww3_dt + !nmth +endif ! WIM + +! ------------------------------------------------------------------------------ - call ice_step + timeLoop: do + +! Noah Day WIM, ---------------------------------------------------------------- + !!!!! LUKE'S WAVE STUFF !!!!! +if (WIM.eq.1) then + ! Update the dates. + + nww3= nww3+nww3_dt + if (WAVE_METH.eq.1) then + !write(nu_diag,*) 'LB: isteep,nww3,N_tm,nmth, mday=', istep, nww3, N_tm, nmth, mday + if (nww3.le.N_tm) then + ww3_swh(:,:) = ww3_swh_full(:,:,nww3) + ww3_fp(:,:) = ww3_fp_full(:,:,nww3) + ww3_dir(:,:) = ww3_dir_full(:,:,nww3) + else + ww3_swh(:,:) = ww3_swh_full(:,:,N_tm) + ww3_fp(:,:) = ww3_fp_full(:,:,N_tm) + ww3_dir(:,:) = ww3_dir_full(:,:,N_tm) + deallocate(ww3_swh_full) + deallocate(ww3_fp_full) + deallocate(ww3_dir_full) + deallocate(ww3_lat) + deallocate(ww3_lon) + deallocate(ww3_tm) + call sub_WW3_dataread(nmth,N_tm,N_lat,N_lon,nyr,mday) + endif + endif + + + ! SETTING CONSTANT VALUES 21/02/23 + if (WAVE_METH.eq.0) then + !allocate(ww3_swh(nx_global,ny_global)) + !allocate(ww3_fp(nx_global,ny_global)) + !allocate(ww3_dir(nx_global,ny_global)) + !ww3_swh(:,:) = 3.0 + !ww3_fp(:,:) = 10.0 + !ww3_dir(:,:) = 0.0 + !write(nu_diag,*) 'Constant wave data set ' + endif + !nmth = nmth+1 + !print*, ' month=', nmth + !if (WAVE_METH.eq.1) call sub_WW3_dataread(nmth,N_tm,N_lat,N_lon) + + !nww3=1-nww3_dt + if (nmth.ne.mmonth) then + nww3 = 24*(mday-1) ! if month changes reset to 1 + nmth = mmonth + write(nu_diag,*) 'New month: isteep,nww3,N_tm,nmth=', istep, nww3, N_tm, nmth + endif + +endif ! WIM + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! ------------------------------------------------------------------------------ + +>>>>>>> Stashed changes +#endif + call ice_step ! tcraig, use advance_timestep now ! istep = istep + 1 ! update time step counters ! istep1 = istep1 + 1 @@ -126,6 +230,195 @@ subroutine CICE_Run call ice_timer_stop(timer_step) ! end timestepping loop timer end subroutine CICE_Run +<<<<<<< Updated upstream +======= +!======================================================================= +!BOP +!! Noah Day WIM +! !ROUTINE: sub_WW3_dataread +! +! !DESCRIPTION: +! +! reads Elodie's WW3 data from netcdf file +! +! !REVISION HISTORY: +! +! authors Luke Bennetts +! +! !INTERFACE: +! + subroutine sub_WW3_dataread (mth,N_tm,N_lat,N_lon,yr,dy) +! +! !USES: +!use ice_read_write, only: ice_open, ice_read, & +! ice_get_ncvarsize, ice_read_vec_nc, & +! ice_open_nc, ice_read_nc, ice_close_nc +use ice_forcing, only: check +use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global +use ice_blocks, only: nx_block, ny_block +! +! !INPUT/OUTPUT PARAMETERS: +! + integer (kind=int_kind), intent(in) :: & + mth,yr,dy ! month 1-12 ! Noah Day 6/3/2022 adding yr, and record number + + integer (kind=int_kind), intent(out) :: & + N_tm, N_lat, N_lon + + ! WW3 variables + integer :: ncid, varid, dimid, numDims, status + integer, dimension(nf90_max_var_dims) :: rhDimIds + + character (char_len_long) :: & ! input data file names + ww3_file, & + varname, & + char_yr, & + tmpname + + integer (kind=int_kind) :: & + fid ! file id for netCDF file + + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: & + work ! output array (real, 8-byte) + +if (WIM.eq.1) then + ! Convert year to character + ! ND: Find what year the model is in. + ! Daily file: + write(tmpname,'(a,a,i4.4,a,a,i4.4,i2.2,i2.2,a)') trim(waveicedatadir), '/', yr, '/', trim(fname_ww3), yr, mth, dy, '.nc' + ! Monthly file: + !write(tmpname,'(a,a,i4.4,a,a,i4.4,i2.2,a)') trim(waveicedatadir), '/', yr, '/', trim(fname_ww3), yr, mth, '.nc' ! Monthly file + + if (WAVE_METH.eq.1) then + if (my_task == master_task) then + write(nu_diag,*) ' sub_WW3_dataread WAVE_METH=1', mth + write(nu_diag,*) ' sub_WW3_dataread file: ', tmpname + !write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.dat' + + endif + ! ND: Read in the WW3 file according to which month the model is in. This is really poorly written but it works. + !call check( nf90_open(trim('/Users/a1724548/GitHub/cice-dirs/input/CICE_data/forcing/gx1/CAWCR/MONTHLY/2005/ww3_200501.nc'), NF90_NOWRITE, ncid )) + if (my_task == master_task) then + !write(nu_diag,*) ' sub_WW3_dataread file: ', trim(waveicedatadir) // '/'// trim(char_yr) //'/'// trim(fname_ww3) // trim(char_yr) // '0101.nc' + write(nu_diag,*) ' sub_WW3_dataread tmpname: ', tmpname + endif + call check( nf90_open(tmpname, NF90_NOWRITE, ncid) ) + + + if (my_task == master_task) then + write(nu_diag,*) '1st check done' + endif + varname = 'LAT'! ND: commenting out 'TLAT' + !varname = 'TLAT' + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + if (my_task == master_task) then + write (nu_diag,*) ' inq is done' + endif + call check( nf90_inquire_variable(ncid, varid, ndims = numDims) ) + call check( nf90_inquire_variable(ncid, varid, dimids = rhDimIds(:numDims)) ) + call check( nf90_inquire_dimension(ncid, rhDimIds(1), len = N_lat) ) + if (my_task == master_task) then + write(nu_diag,*) 'ND: N_lat is:', N_lat + endif + !N_lat = 384 ! grid resolution + N_lat = ny_global !1080 + allocate(ww3_lat(1,N_lat)) ! noah day this used to be (N_lat,1) + + varname = 'LON'! ND: commenting out 'TLON' + !varname = 'TLON' + call check( nf90_get_var(ncid, varid, ww3_lat) ) + if (my_task == master_task) then + write(nu_diag,*) 'Number of dimensions', numDims + write(nu_diag,*) 'lat done', N_lat + write(nu_diag,*) 'min ww3_lat ', minval(ww3_lat) + write(nu_diag,*) 'max ww3_lat ', maxval(ww3_lat) + endif + !write(nu_diag,*) ' ww3_lat(1,N_lat) : ', ww3_lat(1,N_lat) + !write(nu_diag,*) ' ww3_lat(N_lat,1) : ', ww3_lat(N_lat,1) + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + call check( nf90_inquire_variable(ncid, varid, ndims = numDims) ) + call check( nf90_inquire_variable(ncid, varid, dimids = rhDimIds(:numDims)) ) + call check( nf90_inquire_dimension(ncid, rhDimIds(1), len = N_lon) ) + allocate(ww3_lon(N_lon,1)) + if (my_task == master_task) then + write(nu_diag,*) 'lon done', N_lon + endif + varname = 'time' + call check( nf90_get_var(ncid, varid, ww3_lon) ) + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + call check( nf90_inquire_variable(ncid, varid, ndims = numDims) ) + call check( nf90_inquire_variable(ncid, varid, dimids = rhDimIds(:numDims)) ) + call check( nf90_inquire_dimension(ncid, rhDimIds(1), len = N_tm) ) + allocate(ww3_tm(N_tm,1)) + if (my_task == master_task) then + write(nu_diag,*) 'min ww3_lon ', minval(ww3_lon) + write(nu_diag,*) 'max ww3_lon ', maxval(ww3_lon) + write(nu_diag,*) 'time done', N_tm + endif + varname = 'hs' + if (my_task == master_task) then + write(nu_diag,*) 'HS ' + endif + call check( nf90_get_var(ncid, varid, ww3_tm) ) + if (my_task == master_task) then + write(nu_diag,*) 'ww3_tm ', SHAPE(ww3_tm) + endif + allocate(ww3_swh_full(N_lon,N_lat,N_tm)) + if (my_task == master_task) then + write(nu_diag,*) 'allocated', SHAPE(ww3_swh_full) + endif + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + if (my_task == master_task) then + write(nu_diag,*) 'varname: ', trim(varname) + endif + status=nf90_get_var(ncid, varid, ww3_swh_full)!call check( nf90_get_var(ncid, varid, ww3_swh_full) ) + ww3_swh_full = ww3_swh_full!c2*eps3*ww3_swh_full + + !recnum=1 + !call ice_open_nc (trim(waveicedatadir) // '/'// trim(char_yr) //'/'// trim(fname_ww3) // trim(char_yr) // '01.nc', fid) + !call ice_read_nc(ncid,recnum,varname,work(:,:,:)) + !call ice_close_nc(fid) + !ww3_swh_full = c0 + + if (my_task == master_task) then + write(nu_diag,*) 'swh done', SHAPE(ww3_swh_full) + write(nu_diag,*) 'min swh ', minval(ww3_swh_full) + write(nu_diag,*) 'max swh ', maxval(ww3_swh_full) + !write(nu_diag,*) 'swh done', SHAPE(work) + !write(nu_diag,*) 'min swh ', minval(work) + !write(nu_diag,*) 'max swh ', maxval(work) + endif + varname = 'fp' + allocate(ww3_fp_full(N_lon,N_lat,N_tm)) + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + call check( nf90_get_var(ncid, varid, ww3_fp_full) ) + !ww3_fp_full = eps3*ww3_fp_full + ww3_fp_full = ww3_fp_full + if (my_task == master_task) then + write(nu_diag,*) 'fp done', SHAPE(ww3_fp_full) + write(nu_diag,*) 'min fp ', minval(ww3_fp_full) + write(nu_diag,*) 'max fp ', maxval(ww3_fp_full) + endif + varname = 'dir' ! True north degrees, wind direction + allocate(ww3_dir_full(N_lon,N_lat,N_tm)) + call check( nf90_inq_varid(ncid, trim(varname), varid) ) + call check( nf90_get_var(ncid, varid, ww3_dir_full) ) + ww3_dir_full = ww3_dir_full!eps1*ww3_dir_full ! in degrees at this point + if (my_task == master_task) then + write(nu_diag,*) 'dirs done' + write(nu_diag,*) 'min dir ', minval(ww3_dir_full) + write(nu_diag,*) 'max dir ', maxval(ww3_dir_full) + endif + call check( nf90_close(ncid)) + if (my_task == master_task) then + write(nu_diag,*) ' -> WAVE DATA LOADED, N_tm=',N_tm + endif + endif ! WAVE_METH +endif ! WIM + end subroutine sub_WW3_dataread + +!======================================================================= +>>>>>>> Stashed changes !======================================================================= ! diff --git a/configuration/scripts/cice.batch.csh b/configuration/scripts/cice.batch.csh old mode 100755 new mode 100644 index 3e983cd8e..79bf381ae --- a/configuration/scripts/cice.batch.csh +++ b/configuration/scripts/cice.batch.csh @@ -315,6 +315,10 @@ cat >> ${jobfile} << EOFB # nothing to do EOFB +else if (${ICE_MACHINE} =~ nd0349*) then +cat >> ${jobfile} << EOFB +# nothing to do +EOFB else echo "${0} ERROR: ${ICE_MACHINE} unknown" @@ -322,3 +326,4 @@ else endif exit 0 + diff --git a/configuration/scripts/cice.build b/configuration/scripts/cice.build old mode 100755 new mode 100644 diff --git a/configuration/scripts/cice.launch.csh b/configuration/scripts/cice.launch.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/cice.run.setup.csh b/configuration/scripts/cice.run.setup.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/cice.settings b/configuration/scripts/cice.settings old mode 100755 new mode 100644 diff --git a/configuration/scripts/cice.test.setup.csh b/configuration/scripts/cice.test.setup.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/cice_decomp.csh b/configuration/scripts/cice_decomp.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/Macros.nd0349_gnu b/configuration/scripts/machines/Macros.nd0349_gnu new file mode 100644 index 000000000..1bb3e8f5b --- /dev/null +++ b/configuration/scripts/machines/Macros.nd0349_gnu @@ -0,0 +1,47 @@ +#============================================================================== +# Makefile macros for a1724548 - GCC and openmpi compilers +#============================================================================== + +CPP := cpp +CPPDEFS := -DFORTRANUNDERSCORE -DNO_R16 -DHAVE_F2008_CONTIGUOUS -DLINUX -DCPRINTEL ${ICE_CPPDEFS} +CFLAGS := -c -O2 + +FIXEDFLAGS := -132 +FFLAGS := -O2 -ffree-line-length-none -fconvert=big-endian -finit-real=nan +FFLAGS_NOOPT:= -O0 + +ifeq ($(ICE_BLDDEBUG), true) + #FFLAGS += -O0 -g -Wextra -fbacktrace -fbounds-check -ffpe-trap=zero,overflow + FFLAGS += -O0 -g -std=f2008 -fbacktrace -fbounds-check -ffpe-trap=zero,overflow +else + FFLAGS += -O2 +endif + +FC := mpif90 + +MPICC:= + +MPIFC:= mpif90 +LD:= $(FC) + +NETCDF_PATH := $(NETCDF) + +ifeq ($(ICE_IOTYPE), netcdf) + NETCDF_PATH := $(shell nc-config --prefix) + INCLDIR := $(INCLDIR) -I$(NETCDF_PATH)/include + LIB_NETCDF := $(NETCDF_PATH)/lib + LIB_PNETCDF := + SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff +else + SLIBS := +endif + +LIB_MPI := +SCC:=gcc +SFC:= + +#ifeq ($(ICE_THREADED), true) +# LDFLAGS += -fopenmp +# CFLAGS += -fopenmp +# FFLAGS += -fopenmp +#endif diff --git a/configuration/scripts/machines/env.badger_intel b/configuration/scripts/machines/env.badger_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.banting_gnu b/configuration/scripts/machines/env.banting_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.banting_intel b/configuration/scripts/machines/env.banting_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.cesium_intel b/configuration/scripts/machines/env.cesium_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.cheyenne_gnu b/configuration/scripts/machines/env.cheyenne_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.cheyenne_intel b/configuration/scripts/machines/env.cheyenne_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.cheyenne_pgi b/configuration/scripts/machines/env.cheyenne_pgi old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.compy_intel b/configuration/scripts/machines/env.compy_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conda_linux b/configuration/scripts/machines/env.conda_linux old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conda_macos b/configuration/scripts/machines/env.conda_macos old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conrad_cray b/configuration/scripts/machines/env.conrad_cray old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conrad_gnu b/configuration/scripts/machines/env.conrad_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conrad_intel b/configuration/scripts/machines/env.conrad_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.conrad_pgi b/configuration/scripts/machines/env.conrad_pgi old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.cori_intel b/configuration/scripts/machines/env.cori_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.daley_gnu b/configuration/scripts/machines/env.daley_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.daley_intel b/configuration/scripts/machines/env.daley_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.fram_intel b/configuration/scripts/machines/env.fram_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.freya_gnu b/configuration/scripts/machines/env.freya_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.freya_intel b/configuration/scripts/machines/env.freya_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gaffney_gnu b/configuration/scripts/machines/env.gaffney_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gaffney_intel b/configuration/scripts/machines/env.gaffney_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gordon_cray b/configuration/scripts/machines/env.gordon_cray old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gordon_gnu b/configuration/scripts/machines/env.gordon_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gordon_intel b/configuration/scripts/machines/env.gordon_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.gordon_pgi b/configuration/scripts/machines/env.gordon_pgi old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.hera_intel b/configuration/scripts/machines/env.hera_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.high_Sierra_gnu b/configuration/scripts/machines/env.high_Sierra_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.hobart_intel b/configuration/scripts/machines/env.hobart_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.hobart_nag b/configuration/scripts/machines/env.hobart_nag old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.izumi_gnu b/configuration/scripts/machines/env.izumi_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.izumi_intel b/configuration/scripts/machines/env.izumi_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.izumi_nag b/configuration/scripts/machines/env.izumi_nag old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.izumi_pgi b/configuration/scripts/machines/env.izumi_pgi old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.koehr_intel b/configuration/scripts/machines/env.koehr_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.millikan_intel b/configuration/scripts/machines/env.millikan_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.mustang_intel18 b/configuration/scripts/machines/env.mustang_intel18 old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.mustang_intel19 b/configuration/scripts/machines/env.mustang_intel19 old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.mustang_intel20 b/configuration/scripts/machines/env.mustang_intel20 old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.nd0349_gnu b/configuration/scripts/machines/env.nd0349_gnu new file mode 100644 index 000000000..ba023e2c2 --- /dev/null +++ b/configuration/scripts/machines/env.nd0349_gnu @@ -0,0 +1,15 @@ +#!/bin/csh -f + +setenv ICE_MACHINE_MACHNAME nd0349 +setenv ICE_MACHINE_ENVNAME gnu +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR /g/data/ia40/cice-dirs/runs +setenv ICE_MACHINE_INPUTDATA /g/data/ia40/cice-dirs/input +setenv ICE_MACHINE_BASELINE /g/data/ia40/cice-dirs/baseline +setenv ICE_MACHINE_SUBMIT " " +setenv ICE_MACHINE_TPNODE 4 +setenv ICE_MACHINE_ACCT P0000000 +setenv ICE_MACHINE_BLDTHRDS 1 +setenv ICE_MACHINE_QSTAT " " +setenv ICE_MACHINE_QUIETMODE false +setenv ICE_MACHINE_QUEUE " " diff --git a/configuration/scripts/machines/env.onyx_cray b/configuration/scripts/machines/env.onyx_cray old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.onyx_gnu b/configuration/scripts/machines/env.onyx_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.onyx_intel b/configuration/scripts/machines/env.onyx_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.orion_intel b/configuration/scripts/machines/env.orion_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.phase3_intel b/configuration/scripts/machines/env.phase3_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.testmachine_intel b/configuration/scripts/machines/env.testmachine_intel old mode 100755 new mode 100644 diff --git a/configuration/scripts/machines/env.travisCI_gnu b/configuration/scripts/machines/env.travisCI_gnu old mode 100755 new mode 100644 diff --git a/configuration/scripts/makdep.c b/configuration/scripts/makdep.c old mode 100755 new mode 100644 diff --git a/configuration/scripts/options/set_nml.wimom21deg b/configuration/scripts/options/set_nml.wimom21deg new file mode 100644 index 000000000..793474f8e --- /dev/null +++ b/configuration/scripts/options/set_nml.wimom21deg @@ -0,0 +1,31 @@ +tr_fsd = .true. +restart_fsd = .false. +nfsd = 12 +wave_spec_type = 'profile' +nfreq = 31 +ice_ic = '/Users/noahday/GitHub/cice-dirs/input/CICE_data/ic/access-om2_1deg/iced.2019-01-01-00000.nc' +nprocs = 1 +nx_global = 360 +ny_global = 300 +block_size_x = 360 +block_size_y = 20 +max_blocks = 32 +processor_shape = 'slenderX1' +ssh_stress = 'geostrophic' +grid_format = 'nc' +grid_type = 'tripole' +grid_file = '/Users/noahday/GitHub/cice-dirs/input/CICE_data/grid/om2_1deg/icegrid_nonc.nc' +kmt_file = '/Users/noahday/GitHub/cice-dirs/input/CICE_data/grid/om2_1deg/ocean_grid_10.nc' +tfrz_option = 'mushy' +atm_data_type = 'JRA55_gx1' +fyear_init = 2017 +ycycle = 1 +atm_data_format = 'nc' +atm_data_dir = '/Users/noahday/GitHub/cice-dirs/input/CICE_data/forcing/access-om2_1deg/JRA55-do/1-4-0' +ocn_data_format = 'bin' +ocn_data_dir = '/Users/noahday/GitHub/cice-dirs/input/CICE_data/forcing/access-om2_1deg/ocean' +oceanmixed_file = 'ocean_output2017.nc' +ns_boundary_type = 'tripole' +ocn_data_format = 'nc' +nilyr = 4 +ocn_data_format = 'nc' \ No newline at end of file diff --git a/configuration/scripts/parse_namelist.sh b/configuration/scripts/parse_namelist.sh old mode 100755 new mode 100644 diff --git a/configuration/scripts/parse_namelist_from_env.sh b/configuration/scripts/parse_namelist_from_env.sh old mode 100755 new mode 100644 diff --git a/configuration/scripts/parse_settings.sh b/configuration/scripts/parse_settings.sh old mode 100755 new mode 100644 diff --git a/configuration/scripts/set_version_number.csh b/configuration/scripts/set_version_number.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/setup_run_dirs.csh b/configuration/scripts/setup_run_dirs.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/QC/cice.t-test.py b/configuration/scripts/tests/QC/cice.t-test.py old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/QC/compare_qc_cases.csh b/configuration/scripts/tests/QC/compare_qc_cases.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/QC/gen_qc_cases.csh b/configuration/scripts/tests/QC/gen_qc_cases.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/cice_test_codecov.csh b/configuration/scripts/tests/cice_test_codecov.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/comparebfb.csh b/configuration/scripts/tests/comparebfb.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/comparelog.csh b/configuration/scripts/tests/comparelog.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/lcov_modify_source.sh b/configuration/scripts/tests/lcov_modify_source.sh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/poll_queue.csh b/configuration/scripts/tests/poll_queue.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/tests/report_results.csh b/configuration/scripts/tests/report_results.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/timeseries.csh b/configuration/scripts/timeseries.csh old mode 100755 new mode 100644 diff --git a/configuration/scripts/timeseries.py b/configuration/scripts/timeseries.py old mode 100755 new mode 100644 diff --git a/configuration/tools/jra55_datasets/interp_jra55_ncdf_bilinear.py b/configuration/tools/jra55_datasets/interp_jra55_ncdf_bilinear.py old mode 100755 new mode 100644 diff --git a/configuration/tools/jra55_datasets/make_forcing.csh b/configuration/tools/jra55_datasets/make_forcing.csh old mode 100755 new mode 100644 diff --git a/doc/source/user_guide/figures/CICE_Bgrid.png b/doc/source/user_guide/figures/CICE_Bgrid.png old mode 100755 new mode 100644 diff --git a/icepack b/icepack index a80472b54..299ebd51b 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit a80472b547aa6d7a85f8ae5e1449273a323e0371 +Subproject commit 299ebd51b7956c32b63b24ebfb3c435524d6a692