diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 27553c01a..687bf8478 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -83,11 +83,13 @@ subroutine eap (dt) use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn use ice_dyn_shared, only: fcor_blk, ndte, dtei, a_min, m_min, & cosw, sinw, denom1, uvel_init, vvel_init, arlx1i, & - evp_prep1, evp_prep2, stepu, evp_finish + evp_prep1, evp_prep2, stepu, evp_finish, & + basal_stress_coeff, l_basalstress use ice_flux, only: rdg_conv, rdg_shear, prs_sig, strairxT, strairyT, & strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & strtltx, strtlty, strocnx, strocny, strintx, strinty, & strocnxT, strocnyT, strax, stray, & + Cbu, hwater, tau_bu, tau_bv, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 @@ -269,7 +271,8 @@ subroutine eap (dt) stress12_1(:,:,iblk), stress12_2(:,:,iblk), & stress12_3(:,:,iblk), stress12_4(:,:,iblk), & uvel_init (:,:,iblk), vvel_init (:,:,iblk), & - uvel (:,:,iblk), vvel (:,:,iblk)) + uvel (:,:,iblk), vvel (:,:,iblk), & + Cbu (:,:,iblk)) !----------------------------------------------------------------- ! Initialize structure tensor @@ -387,6 +390,21 @@ subroutine eap (dt) strtmp (:,:,:)) ! call ice_timer_stop(timer_tmp1) ! dynamics + !----------------------------------------------------------------- + ! basal stress calculation (landfast ice) + !----------------------------------------------------------------- + + if (l_basalstress) then + call basal_stress_coeff (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), & + uvel(:,:,iblk), vvel(:,:,iblk), & + Cbu(:,:,iblk)) + endif + + !----------------------------------------------------------------- ! momentum equation !----------------------------------------------------------------- @@ -403,7 +421,8 @@ subroutine eap (dt) strocnx (:,:,iblk), strocny (:,:,iblk), & strintx (:,:,iblk), strinty (:,:,iblk), & uvel_init(:,:,iblk), vvel_init(:,:,iblk),& - uvel (:,:,iblk), vvel (:,:,iblk)) + uvel (:,:,iblk), vvel (:,:,iblk), & + Cbu (:,:,iblk)) ! load velocity into array for boundary updates fld2(:,:,1,iblk) = uvel(:,:,iblk) @@ -453,6 +472,16 @@ subroutine eap (dt) !$OMP END PARALLEL DO enddo ! subcycling + + ! calculate basal stress component for outputs + if ( l_basalstress ) then + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + tau_bu(:,:,iblk) = Cbu(:,:,iblk)*uvel(:,:,iblk) + tau_bv(:,:,iblk) = Cbu(:,:,iblk)*vvel(:,:,iblk) + enddo + !$OMP END PARALLEL DO + endif deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 95ce92e8e..13e7b949c 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -37,7 +37,7 @@ module ice_dyn_evp use ice_kinds_mod use ice_dyn_shared, only: stepu, evp_prep1, evp_prep2, evp_finish, & ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, & - vvel_init + vvel_init, basal_stress_coeff, l_basalstress, Ktens implicit none private @@ -77,6 +77,7 @@ subroutine evp (dt) strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & strtltx, strtlty, strocnx, strocny, strintx, strinty, & strocnxT, strocnyT, strax, stray, & + Cbu, hwater, tau_bu, tau_bv, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 @@ -263,7 +264,8 @@ subroutine evp (dt) stress12_1(:,:,iblk), stress12_2(:,:,iblk), & stress12_3(:,:,iblk), stress12_4(:,:,iblk), & uvel_init (:,:,iblk), vvel_init (:,:,iblk), & - uvel (:,:,iblk), vvel (:,:,iblk)) + uvel (:,:,iblk), vvel (:,:,iblk), & + Cbu (:,:,iblk)) !----------------------------------------------------------------- ! ice strength @@ -347,6 +349,20 @@ subroutine evp (dt) strtmp (:,:,:) ) ! endif ! yield_curve + !----------------------------------------------------------------- + ! basal stress calculation (landfast ice) + !----------------------------------------------------------------- + + if (l_basalstress) then + call basal_stress_coeff (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), & + uvel(:,:,iblk), vvel(:,:,iblk), & + Cbu(:,:,iblk)) + endif + !----------------------------------------------------------------- ! momentum equation !----------------------------------------------------------------- @@ -363,7 +379,8 @@ subroutine evp (dt) strocnx (:,:,iblk), strocny (:,:,iblk), & strintx (:,:,iblk), strinty (:,:,iblk), & uvel_init(:,:,iblk), vvel_init(:,:,iblk),& - uvel (:,:,iblk), vvel (:,:,iblk)) + uvel (:,:,iblk), vvel (:,:,iblk), & + Cbu (:,:,iblk)) ! load velocity into array for boundary updates fld2(:,:,1,iblk) = uvel(:,:,iblk) @@ -388,8 +405,18 @@ subroutine evp (dt) vvel(:,:,iblk) = fld2(:,:,2,iblk) enddo !$OMP END PARALLEL DO - + enddo ! subcycling + + ! calculate basal stress component for outputs + if ( l_basalstress ) then + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + tau_bu(:,:,iblk) = Cbu(:,:,iblk)*uvel(:,:,iblk) + tau_bv(:,:,iblk) = Cbu(:,:,iblk)*vvel(:,:,iblk) + enddo + !$OMP END PARALLEL DO + endif deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) @@ -521,7 +548,7 @@ subroutine stress (nx_block, ny_block, & str ) use ice_constants, only: c0, c4, p027, p055, p111, p166, & - p2, p222, p25, p333, p5, puny + p2, p222, p25, p333, p5, puny, c1 integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -683,24 +710,24 @@ subroutine stress (nx_block, ny_block, & ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune - Deltane)) & + stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune*(c1+Ktens) - Deltane*(c1-Ktens))) & * denom1 - stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw - Deltanw)) & + stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens))) & * denom1 - stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw - Deltasw)) & + stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens))) & * denom1 - stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse - Deltase)) & + stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse*(c1+Ktens) - Deltase*(c1-Ktens))) & * denom1 - stressm_1(i,j) = (stressm_1(i,j) + c0ne*tensionne) * denom1 - stressm_2(i,j) = (stressm_2(i,j) + c0nw*tensionnw) * denom1 - stressm_3(i,j) = (stressm_3(i,j) + c0sw*tensionsw) * denom1 - stressm_4(i,j) = (stressm_4(i,j) + c0se*tensionse) * denom1 - - stress12_1(i,j) = (stress12_1(i,j) + c0ne*shearne*p5) * denom1 - stress12_2(i,j) = (stress12_2(i,j) + c0nw*shearnw*p5) * denom1 - stress12_3(i,j) = (stress12_3(i,j) + c0sw*shearsw*p5) * denom1 - stress12_4(i,j) = (stress12_4(i,j) + c0se*shearse*p5) * denom1 + stressm_1(i,j) = (stressm_1(i,j) + c0ne*tensionne*(c1+Ktens)) * denom1 + stressm_2(i,j) = (stressm_2(i,j) + c0nw*tensionnw*(c1+Ktens)) * denom1 + stressm_3(i,j) = (stressm_3(i,j) + c0sw*tensionsw*(c1+Ktens)) * denom1 + stressm_4(i,j) = (stressm_4(i,j) + c0se*tensionse*(c1+Ktens)) * denom1 + + stress12_1(i,j) = (stress12_1(i,j) + c0ne*shearne*p5*(c1+Ktens)) * denom1 + stress12_2(i,j) = (stress12_2(i,j) + c0nw*shearnw*p5*(c1+Ktens)) * denom1 + stress12_3(i,j) = (stress12_3(i,j) + c0sw*shearsw*p5*(c1+Ktens)) * denom1 + stress12_4(i,j) = (stress12_4(i,j) + c0se*shearse*p5*(c1+Ktens)) * denom1 !----------------------------------------------------------------- ! Eliminate underflows. diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 81a5b08de..434ac44fa 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -18,7 +18,8 @@ module ice_dyn_shared implicit none private public :: init_evp, set_evp_parameters, stepu, principal_stress, & - evp_prep1, evp_prep2, evp_finish + evp_prep1, evp_prep2, evp_finish, basal_stress_coeff, & + read_basalstress_bathy save ! namelist parameters @@ -45,6 +46,7 @@ module ice_dyn_shared real (kind=dbl_kind), public :: & revp , & ! 0 for classic EVP, 1 for revised EVP + e_ratio , & ! e = EVP ellipse aspect ratio ecci , & ! 1/e^2 dtei , & ! 1/dte, where dte is subcycling timestep (1/s) dte2T , & ! dte/2T @@ -61,6 +63,13 @@ module ice_dyn_shared uvel_init, & ! x-component of velocity (m/s), beginning of timestep vvel_init ! y-component of velocity (m/s), beginning of timestep + logical (kind=log_kind), public :: & + l_basalstress ! if true, basal stress for landfast on + + ! ice isotropic tensile strength parameter + real (kind=dbl_kind), public :: & + Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) + !======================================================================= contains @@ -184,9 +193,9 @@ subroutine set_evp_parameters (dt) dtei = c1/dte ! 1/s ! major/minor axis length ratio, squared - ecc = c4 - ecci = p25 ! 1/ecc - + ecc = e_ratio**2 + ecci = c1/ecc ! 1/ecc + ! constants for stress equation tdamp2 = c2*eyc*dt ! s dte2T = dte/tdamp2 ! ellipse (unitless) @@ -373,7 +382,8 @@ subroutine evp_prep2 (nx_block, ny_block, & stress12_1, stress12_2, & stress12_3, stress12_4, & uvel_init, vvel_init, & - uvel, vvel) + uvel, vvel, & + Cbu) use ice_constants, only: c0, c1, gravit @@ -420,6 +430,7 @@ subroutine evp_prep2 (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(out) :: & + Cbu, & ! coefficient for basal stress uvel_init,& ! x-component of velocity (m/s), beginning of time step vvel_init,& ! y-component of velocity (m/s), beginning of time step umassdti, & ! mass of U-cell/dt (kg/m^2 s) @@ -462,6 +473,7 @@ subroutine evp_prep2 (nx_block, ny_block, & forcex (i,j) = c0 forcey (i,j) = c0 umassdti (i,j) = c0 + Cbu (i,j) = c0 if (revp==1) then ! revised evp stressp_1 (i,j) = c0 @@ -600,7 +612,8 @@ subroutine stepu (nx_block, ny_block, & strocnx, strocny, & strintx, strinty, & uvel_init, vvel_init,& - uvel, vvel) + uvel, vvel, & + Cbu) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -612,6 +625,7 @@ subroutine stepu (nx_block, ny_block, & indxuj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + Cbu, & ! coefficient for basal stress uvel_init,& ! x-component of velocity (m/s), beginning of timestep vvel_init,& ! y-component of velocity (m/s), beginning of timestep aiu , & ! ice fraction on u-grid @@ -675,7 +689,7 @@ subroutine stepu (nx_block, ny_block, & tauy = vrel*watery(i,j) ! ocn stress term ! revp = 0 for classic evp, 1 for revised evp - cca = (brlx + revp)*umassdti(i,j) + vrel * cosw ! kg/m^2 s + cca = (brlx + revp)*umassdti(i,j) + vrel * cosw + Cbu(i,j) ! kg/m^2 s ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s ab2 = cca**2 + ccb**2 @@ -802,6 +816,148 @@ subroutine evp_finish (nx_block, ny_block, & end subroutine evp_finish +!======================================================================= +! Computes basal stress Cb coefficients (landfast ice) +! +! Lemieux, J. F., B. Tremblay, F. Dupont, M. Plante, G. Smith, D. Dumont (2015). +! A basal stress parameterization form modeling landfast ice. J. Geophys. Res. +! Oceans, 120, 3157-3173. +! +! author: Philippe Blain, CMC (coop summer 2015) +! + subroutine basal_stress_coeff (nx_block, ny_block, icellu, & + indxui, indxuj, & + vice, aice, & + hwater, & + uold, vold, & + Cbu) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icellu ! no. of cells where icetmask = 1 + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxui , & ! compressed index in i-direction + indxuj ! compressed index in j-direction + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + aice , & ! concentration of ice at tracer location + vice , & ! volume per unit area of ice at tracer location + hwater , & ! water depth at tracer location + uold , & ! u component of ice speed at previous iteration + vold ! v component of ice speed at previous iteration + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & + Cbu ! coefficient for basal stress + +! +!EOP + + real (kind=dbl_kind) :: & + au, & ! concentration of ice at u location + hu, & ! volume per unit area of ice at u location (mean thickness) + hwu, & ! water depth at u location + hcu, & ! critical thickness at u location + k1 = 8.0_dbl_kind , & ! first free parameter for landfast parametrization + k2 = 15.0_dbl_kind, & ! second free parameter (Nm^-3) for landfast parametrization + u0 = 5e-5_dbl_kind, & ! residual velocity (m/s) + CC = 20.0_dbl_kind ! CC=Cb factor in Lemieux et al 2015 + + integer (kind=int_kind) :: & + i, j, ij + + do ij = 1, icellu + i = indxui(ij) + j = indxuj(ij) + + ! convert quantities to u-location + au = max(aice(i,j),aice(i+1,j),aice(i,j+1),aice(i+1,j+1)) + hwu = min(hwater(i,j),hwater(i+1,j),hwater(i,j+1),hwater(i+1,j+1)) + hu = max(vice(i,j),vice(i+1,j),vice(i,j+1),vice(i+1,j+1)) + + ! calculate basal stress factor + ! 1- calculate critical thickness + hcu = au * hwu / k1 + + ! 2- calculate stress factor + if (au > p01 .and. hu > hcu ) then + ! endif + Cbu(i,j) = ( k2 / (sqrt(uold(i,j)**2 + vold(i,j)**2) + u0) ) & + * (hu - hcu) * exp(-CC * (1 - au)) + endif + + enddo ! ij + + end subroutine basal_stress_coeff + +!======================================================================= + +! Read bathymetry data for basal stress calculation (grounding scheme for +! landfast ice) in CICE stand-alone mode. When CICE is in coupled mode +! (e.g. CICE-NEMO), hwater should be uptated at each time level so that +! it varies with ocean dynamics. +! +! author: Fred Dupont, CMC + + subroutine read_basalstress_bathy + + ! use module + use ice_blocks, only: block, get_block, nx_block, ny_block + use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn + use ice_domain_size, only: max_blocks + use ice_flux, only: hwater + use ice_read_write + use ice_fileunits, only: nu_diag + use ice_communicate, only: my_task, master_task + use ice_constants, only: field_loc_center, field_type_scalar + + + ! local variables + integer (kind=int_kind) :: & + i, j, & ! index inside block + iblk, & ! block index + fid_init ! file id for netCDF init file + + character (char_len_long) :: & ! input data file names + init_file, & + fieldname + + logical (kind=log_kind) :: diag=.true. + + init_file='bathymetry.nc' + + if (my_task == master_task) then + + write (nu_diag,*) ' ' + write (nu_diag,*) 'Initial ice file: ', trim(init_file) + write (*,*) 'Initial ice file: ', trim(init_file) + call flush(nu_diag) + + endif + + call ice_open_nc(init_file,fid_init) + + fieldname='Bathymetry' + + if (my_task == master_task) then + write(nu_diag,*) 'reading ',TRIM(fieldname) + write(*,*) 'reading ',TRIM(fieldname) + call flush(nu_diag) + endif + call ice_read_nc(fid_init,1,fieldname,hwater,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + call ice_close_nc(fid_init) + + if (my_task == master_task) then + write(nu_diag,*) 'closing file ',TRIM(init_file) + call flush(nu_diag) + endif + + end subroutine read_basalstress_bathy + !======================================================================= ! Computes principal stresses for comparison with the theoretical diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index 53f1acd34..b3c5ff7e2 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -41,6 +41,7 @@ module ice_flux vocn , & ! ocean current, y-direction (m/s) ss_tltx , & ! sea surface slope, x-direction (m/m) ss_tlty , & ! sea surface slope, y-direction + hwater , & ! water depth for basal stress calc (landfast ice) ! out to atmosphere strairxT, & ! stress on ice by air, x-direction @@ -56,6 +57,8 @@ module ice_flux real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: & sig1 , & ! principal stress component sig2 , & ! principal stress component + tau_bu , & ! basal stress (x) (N/m^2) + tau_bv , & ! basal stress (y) (N/m^2) strairx , & ! stress on ice by air, x-direction strairy , & ! stress on ice by air, y-direction strocnx , & ! ice-ocean stress, x-direction @@ -103,7 +106,8 @@ module ice_flux real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: & prs_sig , & ! replacement pressure, for stress calc - fm ! Coriolis param. * mass in U-cell (kg/s) + fm , & ! Coriolis param. * mass in U-cell (kg/s) + Cbu ! coefficient for basal stress (landfast ice) !----------------------------------------------------------------- ! Thermodynamic component @@ -426,6 +430,7 @@ subroutine init_coupler_flux vocn (:,:,:) = c0 frzmlt(:,:,:) = c0 ! freezing/melting potential (W/m^2) sss (:,:,:) = 34.0_dbl_kind ! sea surface salinity (ppt) + hwater(:,:,:) = 10000.0_dbl_kind! water depth for basal stress calc (landfast ice) do iblk = 1, size(Tf,3) do j = 1, size(Tf,2) @@ -653,6 +658,8 @@ subroutine init_history_dyn sig1 (:,:,:) = c0 sig2 (:,:,:) = c0 + tau_bu (:,:,:) = c0 + tau_bv (:,:,:) = c0 strocnx (:,:,:) = c0 strocny (:,:,:) = c0 strairx (:,:,:) = c0 diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 43e14d2c8..e8450881a 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -74,7 +74,8 @@ subroutine input_data sss_data_type, sst_data_type, ocn_data_dir, & oceanmixed_file, restore_sst, trestore use ice_grid, only: grid_file, gridcpl_file, kmt_file, grid_type, grid_format - use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve + use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, & + l_basalstress, Ktens, e_ratio use ice_transport_driver, only: advection use icepack_intfc_tracers, only: tr_iage, tr_FY, tr_lvl, tr_pond, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo, & @@ -133,7 +134,7 @@ subroutine input_data kdyn, ndte, revised_evp, yield_curve, & advection, & kstrength, krdg_partic, krdg_redist, mu_rdg, & - Cf + e_ratio, Ktens, Cf, l_basalstress namelist /shortwave_nml/ & shortwave, albedo_type, & @@ -227,6 +228,9 @@ subroutine input_data krdg_redist = 1 ! 1 = new redistribution, 0 = Hibler 80 mu_rdg = 3 ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5) Cf = 17.0_dbl_kind ! ratio of ridging work to PE change in ridging + l_basalstress= .false. ! if true, basal stress for landfast is on + Ktens = 0.0_dbl_kind ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) + e_ratio = 2.0_dbl_kind ! EVP ellipse aspect ratio advection = 'remap' ! incremental remapping transport scheme shortwave = 'default' ! 'default' or 'dEdd' (delta-Eddington) albedo_type = 'default'! or 'constant' @@ -685,6 +689,9 @@ subroutine input_data call broadcast_scalar(krdg_redist, master_task) call broadcast_scalar(mu_rdg, master_task) call broadcast_scalar(Cf, master_task) + call broadcast_scalar(l_basalstress, master_task) + call broadcast_scalar(Ktens, master_task) + call broadcast_scalar(e_ratio, master_task) call broadcast_scalar(advection, master_task) call broadcast_scalar(shortwave, master_task) call broadcast_scalar(albedo_type, master_task) @@ -857,6 +864,9 @@ subroutine input_data write(nu_diag,1000) ' mu_rdg = ', mu_rdg if (kstrength == 1) & write(nu_diag,1000) ' Cf = ', Cf + write(nu_diag,1010) ' l_basalstress = ', l_basalstress + write(nu_diag,1005) ' Ktens = ', Ktens + write(nu_diag,1005) ' e_ratio = ', e_ratio write(nu_diag,1030) ' advection = ', & trim(advection) write(nu_diag,1030) ' shortwave = ', & diff --git a/cicecore/drivers/cice/CICE_InitMod.F90 b/cicecore/drivers/cice/CICE_InitMod.F90 index 37320b3df..514826921 100644 --- a/cicecore/drivers/cice/CICE_InitMod.F90 +++ b/cicecore/drivers/cice/CICE_InitMod.F90 @@ -62,7 +62,8 @@ subroutine cice_init use ice_domain, only: init_domain_blocks use ice_domain_size, only: ncat use ice_dyn_eap, only: init_eap - use ice_dyn_shared, only: kdyn, init_evp + use ice_dyn_shared, only: kdyn, init_evp, & + l_basalstress, read_basalstress_bathy use ice_fileunits, only: init_fileunits, nu_diag use ice_flux, only: init_coupler_flux, init_history_therm, & init_history_dyn, init_flux_atm, init_flux_ocn @@ -127,6 +128,8 @@ subroutine cice_init call calendar(time) ! determine the initial date call init_forcing_ocn(dt) ! initialize sss and sst from data + if (l_basalstress) & + call read_basalstress_bathy ! read bathy for basalstress calc (standalone mode) call init_state ! initialize the ice state call init_transport ! initialize horizontal transport call ice_HaloRestore_init ! restored boundary conditions diff --git a/cicecore/drivers/hadgem3/CICE_InitMod.F90 b/cicecore/drivers/hadgem3/CICE_InitMod.F90 index a3a3927e2..d6265d05e 100644 --- a/cicecore/drivers/hadgem3/CICE_InitMod.F90 +++ b/cicecore/drivers/hadgem3/CICE_InitMod.F90 @@ -61,7 +61,8 @@ subroutine cice_init use ice_domain, only: init_domain_blocks use ice_domain_size, only: ncat use ice_dyn_eap, only: init_eap - use ice_dyn_shared, only: kdyn, init_evp + use ice_dyn_shared, only: kdyn, init_evp, & + l_basalstress, read_basalstress_bathy use ice_fileunits, only: init_fileunits, nu_diag use ice_flux, only: init_coupler_flux, init_history_therm, & init_history_dyn, init_flux_atm, init_flux_ocn @@ -83,7 +84,7 @@ subroutine cice_init #ifdef popcice use drv_forcing, only: sst_sss #endif - + call init_communicate ! initial setup for message passing call init_fileunits ! unit numbers call input_data ! namelist variables @@ -117,6 +118,8 @@ subroutine cice_init #ifndef CICE_IN_NEMO call init_forcing_ocn(dt) ! initialize sss and sst from data + if (l_basalstress) & + call read_basalstress_bathy ! read bathy for basalstress calculation (standalone mode). #endif call init_state ! initialize the ice state call init_transport ! initialize horizontal transport diff --git a/configuration/data/gx3/bathymetry.nc b/configuration/data/gx3/bathymetry.nc new file mode 100644 index 000000000..7c1295a30 Binary files /dev/null and b/configuration/data/gx3/bathymetry.nc differ diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 167c6a3a0..9e7a93876 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -100,6 +100,9 @@ krdg_redist = 1 mu_rdg = 3 Cf = 17. + e_ratio = 2. + Ktens = 0. + l_basalstress = .false. / &shortwave_nml diff --git a/configuration/scripts/machines/Macros.fram b/configuration/scripts/machines/Macros.fram new file mode 100755 index 000000000..5181a9e18 --- /dev/null +++ b/configuration/scripts/machines/Macros.fram @@ -0,0 +1,70 @@ +#============================================================================== +# Makefile macros for "fram" +#============================================================================== +# For use with intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE -DNO_R16 -DHAVE_F2008_CONTIGUOUS -DLINUX -DCPRINTEL ${CICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise +#-xHost + +FIXEDFLAGS := -132 +FREEFLAGS := -FR +FFLAGS := -O2 -fp-model precise -convert big_endian -assume byterecl -ftz -traceback +#-xHost +FFLAGS_NOOPT:= -O0 + +ifeq ($(CICE_COMMDIR), mpi) + FC := s.f90 -mpi +else + FC := s.f90 +endif + +MPICC:= s.cc -mpi + +MPIFC:= s.f90 -mpi +LD:= $(MPIFC) + +NETCDF_PATH := $(NETCDF) + +PIO_CONFIG_OPTS:= --enable-filesystem-hints=gpfs + +#PNETCDF_PATH := $(PNETCDF) +#PNETCDF_PATH := /glade/u/apps/ch/opt/pio/2.2/mpt/2.15f/intel/17.0.1/lib + +INCLDIR := $(INCLDIR) + +LIB_NETCDF := $(NETCDF_PATH)/lib +LIB_PNETCDF := $(PNETCDF_PATH)/lib +LIB_MPI := $(IMPILIBDIR) + +#SLIBS := -L$(LIB_NETCDF) -lnetcdff -lnetcdf -L$(LIB_PNETCDF) -lpnetcdf -lgptl +SLIBS := -L$(LIB_NETCDF) -lnetcdff -lnetcdf + +SCC:= s.cc + +SFC:= s.f90 + +ifeq ($(compile_threaded), true) + LDFLAGS += -openmp + CFLAGS += -openmp + FFLAGS += -openmp +endif + +ifeq ($(DITTO), yes) + CPPDEFS := $(CPPDEFS) -DREPRODUCIBLE +endif + +### if using parallel I/O, load all 3 libraries. PIO must be first! +ifeq ($(IO_TYPE), pio) + PIO_PATH:=/glade/u/apps/ch/opt/pio/2.2/mpt/2.15f/intel/17.0.1/lib + INCLDIR += -I/glade/u/apps/ch/opt/pio/2.2/mpt/2.15f/intel/17.0.1/include + SLIBS := $(SLIBS) -L$(PIO_PATH) -lpiof + + CPPDEFS := $(CPPDEFS) -Dncdf +endif + +ifeq ($(IO_TYPE), netcdf) + CPPDEFS := $(CPPDEFS) -Dncdf +endif diff --git a/configuration/scripts/machines/env.fram b/configuration/scripts/machines/env.fram new file mode 100755 index 000000000..c546cf0b8 --- /dev/null +++ b/configuration/scripts/machines/env.fram @@ -0,0 +1,20 @@ +#!/bin/csh -f + +#. ssmuse-sh -d /fs/ssm/eccc/mrd/rpn/OCEAN/cncpt-3.1.2 +#source NEMO_compiler.ksh + +setenv CICE_MACHINE_ENVNAME fram +setenv CICE_MACHINE_WKDIR /users/dor/armn/jfl/local1/CICE6/tests/CICE_RUNS +setenv CICE_MACHINE_INPUTDATA /fs/cetus/fs1/cmd/e/afsg/fdu/models/fdm/CICE/code/cice_6.0/CICE/configuration/data/gx3Ncar +setenv CICE_MACHINE_BASELINE /users/dor/armn/jfl/local1/CICE6/tests/CICE_BASELINE +setenv CICE_MACHINE_SUBMIT "qsub" +setenv CICE_MACHINE_TPNODE 36 +setenv CICE_MACHINE_ACCT P0000000 +setenv CICE_MACHINE_BLDTHRDS 1 + +if (-e ~/.cice_proj) then + set account_name = `head -1 ~/.cice_proj` + setenv CICE_ACCT ${account_name} +endif + +echo "je suis dans env" \ No newline at end of file