diff --git a/apps/arctic-20km/arctic20km.py b/apps/arctic-20km/arctic20km.py index 07ba4f7..508e019 100644 --- a/apps/arctic-20km/arctic20km.py +++ b/apps/arctic-20km/arctic20km.py @@ -21,8 +21,8 @@ # Choose a predefined ROMS-application: app='arctic-20km' # Arctic-20km -start_date = datetime(1997,01,16,12) -end_date = datetime(1997,01,17,12) +start_date = datetime(1997,01,17,00) +end_date = datetime(1997,02,16,00) a20params=Params(app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=icecpu,restart=False) diff --git a/apps/arctic-20km/include/arctic-20km.h b/apps/arctic-20km/include/arctic-20km.h index f873177..734c82e 100644 --- a/apps/arctic-20km/include/arctic-20km.h +++ b/apps/arctic-20km/include/arctic-20km.h @@ -26,10 +26,10 @@ #define UV_ADV /* turn ON or OFF advection terms */ #define UV_COR /* turn ON or OFF Coriolis term */ -#undef UV_VIS2 /* turn ON or OFF Laplacian horizontal mixing */ -#define UV_VIS4 /* turn ON or OFF biharmonic horizontal mixing */ +#define UV_VIS2 /* turn ON or OFF Laplacian horizontal mixing */ +#undef UV_VIS4 /* turn ON or OFF biharmonic horizontal mixing */ #undef UV_U3ADV_SPLIT /* use 3rd-order upstream split momentum advection */ -#undef UV_U3HADVECTION /* define if 3rd-order upstream horiz. advection */ +#define UV_U3HADVECTION /* define if 3rd-order upstream horiz. advection */ #undef UV_SADVECTION /* turn ON or OFF splines vertical advection */ #undef UV_C4HADVECTION /* define if 4th-order centered horizontal advection */ #define UV_QDRAG /* turn ON or OFF quadratic bottom friction */ @@ -38,7 +38,8 @@ #undef VISC_GRID /* viscosity coefficient scaled by grid size */ #define NONLIN_EOS /* define if using nonlinear equation of state */ #undef WJ_GRADP /* Weighted density Jacobian (Song, 1998) */ -#define DJ_GRADPS /* Splines density Jacobian (Shchepetkin, 2000) */ +#define DJ_GRADPS /* Splines density Jacobian (Shchepetkin, 2000) */ +#undef PJ_GRADPQ4 /* quartic 4 Pressure Jacobian (Shchepetkin,2000) */ #undef DIFF_GRID /* diffusion coefficient scaled by grid size */ #define TS_DIF2 /* turn ON or OFF Laplacian horizontal mixing */ @@ -46,7 +47,7 @@ #undef TS_U3ADV_SPLIT /* use 3rd-order upstream split tracer advection */ #undef TS_U3HADVECTION /* define if 3rd-order upstream horiz. advection */ #define TS_A4HADVECTION /* define if 4th-order Akima horiz. advection */ -#undef TS_C4HADVECTION /* define if 4th-order centered horizontal advection */ +#undef TS_C4HADVECTION /* define if 4th-order centered horizontal advection */ #undef TS_MPDATA /* define if recursive MPDATA 3D advection */ @@ -75,12 +76,15 @@ #undef STATIONS_CGRID /* define if extracting data at native C-grid */ #undef BVF_MIXING /* define if Brunt_Vaisala frequency mixing */ -#define LMD_MIXING /* define if Large et al. (1994) interior closure */ +#define LMD_MIXING /* define if Large et al. (1994) interior closure */ #undef MY25_MIXING /* define if Mellor/Yamada level-2.5 mixing */ #undef GLS_MIXING /* Activate Generic Length-Scale mixing */ #ifdef GLS_MIXING # define N2S2_HORAVG /* Activate horizontal smoothing of buoyancy/shear */ +# undef KANTHA_CLAYSON /* Value for CLS_CMU0 and CLS_C3M vary with choise of stability function */ +# define CANUTO_A +# undef CANUTO_B #endif #ifdef LMD_MIXING # define LMD_BKPP /* use if bottom boundary layer KPP mixing */ diff --git a/apps/build_roms.sh b/apps/build_roms.sh index 6652f95..81fbb04 100755 --- a/apps/build_roms.sh +++ b/apps/build_roms.sh @@ -72,6 +72,7 @@ export USE_OpenMP= export USE_LARGE=on export USE_DEBUG= + export USE_NETCDF4=on export USE_CICE=on diff --git a/apps/common/cice_input_grids/arctic-20km/ice_in.cice5.1.2 b/apps/common/cice_input_grids/arctic-20km/ice_in.cice5.1.2 index e3f83cc..d1f925f 100755 --- a/apps/common/cice_input_grids/arctic-20km/ice_in.cice5.1.2 +++ b/apps/common/cice_input_grids/arctic-20km/ice_in.cice5.1.2 @@ -16,9 +16,9 @@ , restart_dir = '/rundir/restart/' , restart_file = 'iced' , pointer_file = '/rundir/restart/ice.restart_file' - , dumpfreq = 'y' - , dumpfreq_n = 1 - , dump_last = .true. + , dumpfreq = 'd' + , dumpfreq_n = 30 + , dump_last = .false. , bfbflag = .false. , diagfreq = 24 , diag_type = 'stdout' @@ -93,7 +93,7 @@ &dynamics_nml kdyn = 1 , ndte = 120 - , revised_evp = .false. + , revised_evp = .true. , advection = 'remap' , kstrength = 1 , krdg_partic = 1 @@ -152,11 +152,11 @@ &forcing_nml formdrag = .false. , atmbndy = 'default' - , fyear_init = 1997 + , fyear_init = 1993 , ycycle = 100 , atm_data_format = 'nc' - , atm_data_type = 'ecmwf' - , atm_data_dir = './' + , atm_data_type = 'metroms' + , atm_data_dir = './atmcice/' , calc_strair = .true. , highfreq = .false. , natmiter = 5 @@ -164,7 +164,7 @@ , precip_units = 'm_per_12hr' , ustar_min = 0.0005 , fbot_xfer_type = 'constant' - , update_ocn_f = .false. + , update_ocn_f = .true. , l_mpond_fresh = .false. , tfrz_option = 'mushy' , oceanmixed_ice = .false. diff --git a/apps/common/modified_src/cice5.1.2/source/ice_mechred.F90 b/apps/common/modified_src/cice5.1.2/source/ice_mechred.F90 new file mode 100644 index 0000000..70387ec --- /dev/null +++ b/apps/common/modified_src/cice5.1.2/source/ice_mechred.F90 @@ -0,0 +1,2276 @@ +! SVN:$Id: ice_mechred.F90 933 2015-03-16 16:52:38Z eclare $ +!======================================================================= + +! Driver for ice mechanical redistribution (ridging) +! +! See these references: +! +! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength +! in modeling the thickness distribution of Arctic sea ice, +! J. Geophys. Res., 100, 18,611-18,626. +! +! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice +! cover, Mon. Wea. Rev., 108, 1943-1973, 1980. +! +! Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki, 2007: +! Improving ridging schemes for high-resolution sea ice models. +! J. Geophys. Res. 112, C03S91, doi:10.1029/2005JC003355. +! +! Rothrock, D. A., 1975: The energetics of the plastic deformation of +! pack ice by ridging, J. Geophys. Res., 80, 4514-4519. +! +! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, +! 1975: The thickness distribution of sea ice, J. Geophys. Res., +! 80, 4501-4513. +! +! authors: William H. Lipscomb, LANL +! Elizabeth C. Hunke, LANL +! +! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb +! 2004: Block structure added by William Lipscomb +! 2006: New options for participation and redistribution (WHL) +! 2006: Streamlined for efficiency by Elizabeth Hunke +! Converted to free source form (F90) + + module ice_mechred + + use ice_kinds_mod + use ice_constants + use ice_domain_size, only: ncat, max_aero, n_aero, nilyr, nslyr, nblyr + use ice_fileunits, only: nu_diag + use ice_itd, only: hin_max, column_sum, & + column_conservation_check, compute_tracers + + implicit none + save + + private + public :: ice_strength, ridge_ice + +!----------------------------------------------------------------------- +! Ridging parameters +!----------------------------------------------------------------------- + + integer (kind=int_kind), public :: & ! defined in namelist + kstrength , & ! 0 for simple Hibler (1979) formulation + ! 1 for Rothrock (1975) pressure formulation + krdg_partic , & ! 0 for Thorndike et al. (1975) formulation + ! 1 for exponential participation function + krdg_redist ! 0 for Hibler (1980) formulation + ! 1 for exponential redistribution function + + real (kind=dbl_kind), public :: & + mu_rdg, & ! gives e-folding scale of ridged ice (m^.5) (krdg_redist = 1) + Cf ! ratio of ridging work to PE change in ridging (kstrength = 1) + + real (kind=dbl_kind), parameter :: & + Cs = p25 , & ! fraction of shear energy contrbtng to ridging + Cp = p5*gravit*(rhow-rhoi)*rhoi/rhow, & ! proport const for PE + fsnowrdg = p5 , & ! snow fraction that survives in ridging + Gstar = p15 , & ! max value of G(h) that participates + ! (krdg_partic = 0) + astar = p05 , & ! e-folding scale for G(h) participation +!echmod astar = p1 , & ! e-folding scale for G(h) participation + ! (krdg_partic = 1) + maxraft= c1 , & ! max value of hrmin - hi = max thickness + ! of ice that rafts (m) + Hstar = c25 , & ! determines mean thickness of ridged ice (m) + ! (krdg_redist = 0) + ! Flato & Hibler (1995) have Hstar = 100 + Pstar = 2.75e4_dbl_kind, & ! constant in Hibler strength formula + ! (kstrength = 0) + Cstar = c20 ! constant in Hibler strength formula + ! (kstrength = 0) + + logical (kind=log_kind), parameter :: & +! l_conservation_check = .true. ! if true, check conservation + l_conservation_check = .false. ! if true, check conservation + ! (useful for debugging) + +!======================================================================= + + contains + +!======================================================================= + +! Compute changes in the ice thickness distribution due to divergence +! and shear. +! +! author: William H. Lipscomb, LANL + + subroutine ridge_ice (nx_block, ny_block, & + dt, ndtd, & + ntrcr, icells, & + indxi, indxj, & + rdg_conv, rdg_shear, & + aicen, trcrn, & + vicen, vsnon, & + aice0, & + trcr_depend, l_stop, & + istop, jstop, & + dardg1dt, dardg2dt, & + dvirdgdt, opening, & + fpond, & + fresh, fhocn, & + faero_ocn, & + aparticn, krdgn, & + aredistn, vredistn, & + dardg1ndt, dardg2ndt, & + dvirdgndt, & + araftn, vraftn) + + use ice_state, only: nt_qice, nt_qsno, tr_brine, nt_fbri, nt_sice + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells , & ! number of cells with ice present + ndtd , & ! number of dynamics subcycles + ntrcr ! number of tracers in use + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & + rdg_conv, & ! normalized energy dissipation due to convergence (1/s) + rdg_shear ! normalized energy dissipation due to shear (1/s) + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(inout) :: & + aicen , & ! concentration of ice + vicen , & ! volume per unit area of ice (m) + vsnon ! volume per unit area of snow (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), & + intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), dimension (nx_block,ny_block), & + intent(inout) :: & + aice0 ! concentration of open water + + integer (kind=int_kind), dimension(ntrcr), intent(in) :: & + trcr_depend + + logical (kind=log_kind), intent(out) :: & + l_stop ! if true, abort on return + + integer (kind=int_kind), intent(out) :: & + istop, jstop ! indices of grid cell where model aborts + + ! optional history fields + real (kind=dbl_kind), dimension(nx_block,ny_block), & + intent(inout), optional :: & + dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2dt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgdt , & ! rate of ice volume ridged (m/s) + opening , & ! rate of opening due to divergence/shear (1/s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), & + intent(out), optional :: & + dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgndt , & ! rate of ice volume ridged (m/s) + aparticn , & ! participation function + krdgn , & ! mean ridge thickness/thickness of ridging ice + araftn , & ! rafting ice area + vraftn , & ! rafting ice volume + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume + + real (kind=dbl_kind), dimension(nx_block,ny_block,max_aero), & + intent(inout), optional :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + ! local variables + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat) :: & + eicen ! energy of melting for each ice layer (J/m^2) + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat) :: & + esnon , & ! energy of melting for each snow layer (J/m^2) + vbrin, & ! ice volume with defined by brine height (m) + sicen ! Bulk salt in h ice (ppt*m) + + real (kind=dbl_kind), dimension (icells) :: & + asum , & ! sum of ice and open water area + aksum , & ! ratio of area removed to area ridged + msnow_mlt , & ! mass of snow added to ocean (kg m-2) + esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) + mpond , & ! mass of pond added to ocean (kg m-2) + closing_net, & ! net rate at which area is removed (1/s) + ! (ridging ice area - area of new ridges) / dt + divu_adv , & ! divu as implied by transport scheme (1/s) + opning , & ! rate of opening due to divergence/shear + ! opning is a local variable; + ! opening is the history diagnostic variable + ardg1 , & ! fractional area loss by ridging ice + ardg2 , & ! fractional area gain by new ridges + virdg , & ! ice volume ridged + aopen ! area opening due to divergence/shear + + real (kind=dbl_kind), dimension (icells,max_aero) :: & + maero ! aerosol mass added to ocean (kg m-2) + + real (kind=dbl_kind), dimension (icells,0:ncat) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (icells,ncat) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg , & ! mean ridge thickness/thickness of ridging ice + ardg1n , & ! area of ice ridged + ardg2n , & ! area of new ridges + virdgn , & ! ridging ice volume + mraftn ! rafting ice mask + + real (kind=dbl_kind), dimension (icells) :: & + vice_init, vice_final, & ! ice volume summed over categories + vsno_init, vsno_final, & ! snow volume summed over categories + eice_init, eice_final, & ! ice energy summed over layers + vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories + Shi_init , Shi_final , & ! ice bulk salinity summed over categories + esno_init, esno_final ! snow energy summed over layers + + integer (kind=int_kind), parameter :: & + nitermax = 20 ! max number of ridging iterations + + integer (kind=int_kind) :: & + i,j , & ! horizontal indices + n , & ! thickness category index + niter , & ! iteration counter + ij , & ! horizontal index, combines i and j loops + k ! vertical index + + real (kind=dbl_kind) :: & + dti ! 1 / dt + + logical (kind=log_kind) :: & + iterate_ridging ! if true, repeat the ridging + + character (len=char_len) :: & + fieldid ! field identifier + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + + l_stop = .false. + istop = 0 + jstop = 0 + + do ij = 1, icells + msnow_mlt(ij) = c0 + esnow_mlt(ij) = c0 + maero (ij,:) = c0 + mpond (ij) = c0 + ardg1 (ij) = c0 + ardg2 (ij) = c0 + virdg (ij) = c0 + ardg1n (ij,:) = c0 + ardg2n (ij,:) = c0 + virdgn (ij,:) = c0 + mraftn (ij,:) = c0 +! aopen (ij) = c0 + enddo + + !----------------------------------------------------------------- + ! Compute area of ice plus open water before ridging. + !----------------------------------------------------------------- + call asum_ridging (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, aice0, & + asum) + + !----------------------------------------------------------------- + ! Compute the area opening and closing. + !----------------------------------------------------------------- + call ridge_prep (nx_block, ny_block, & + icells, indxi, indxj, & + dt, & + rdg_conv, rdg_shear, & + asum, closing_net, & + divu_adv, opning) + + !----------------------------------------------------------------- + ! Compute initial values of conserved quantities. + !----------------------------------------------------------------- + + if (l_conservation_check) then + + eicen(:,:,:) = c0 + esnon(:,:,:) = c0 + vbrin(:,:,:) = c0 + sicen(:,:,:) = c0 + + do n = 1, ncat + do j = 1, ny_block + do i = 1, nx_block + vbrin(i,j,n) = vicen(i,j,n) + if (tr_brine) vbrin(i,j,n) = trcrn(i,j,nt_fbri,n) * vicen(i,j,n) + enddo + enddo + + do k = 1, nilyr + do j = 1, ny_block + do i = 1, nx_block + eicen(i,j,n) = eicen(i,j,n) + trcrn(i,j,nt_qice+k-1,n) & + * vicen(i,j,n)/real(nilyr,kind=dbl_kind) + enddo + enddo + enddo + do k = 1, nslyr + do j = 1, ny_block + do i = 1, nx_block + esnon(i,j,n) = esnon(i,j,n) + trcrn(i,j,nt_qsno+k-1,n) & + * vsnon(i,j,n)/real(nslyr,kind=dbl_kind) + enddo + enddo + enddo + + + do k = 1, nilyr + do j = 1, ny_block + do i = 1, nx_block + sicen(i,j,n) = sicen(i,j,n) + trcrn(i,j,nt_sice+k-1,n) & + * vicen(i,j,n)/real(nilyr,kind=dbl_kind) + enddo + enddo + enddo + + enddo + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vicen, vice_init) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vsnon, vsno_init) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + eicen, eice_init) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + esnon, esno_init) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vbrin, vbri_init) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + sicen, Shi_init) + + endif + + do niter = 1, nitermax + + !----------------------------------------------------------------- + ! Compute the thickness distribution of ridging ice + ! and various quantities associated with the new ridged ice. + !----------------------------------------------------------------- + + call ridge_itd (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, vicen, & + aice0, & + aksum, apartic, & + hrmin, hrmax, & + hrexp, krdg, & + aparticn, krdgn, mraftn) + + !----------------------------------------------------------------- + ! Redistribute area, volume, and energy. + !----------------------------------------------------------------- + + call ridge_shift (nx_block, ny_block, & + icells, indxi, indxj, & + ntrcr, dt, & + aicen, trcrn, & + vicen, vsnon, & + aice0, trcr_depend, & + aksum, apartic, & + hrmin, hrmax, & + hrexp, krdg, & + closing_net, opning, & + ardg1, ardg2, & + virdg, aopen, & + ardg1n, ardg2n, & + virdgn, & + msnow_mlt, esnow_mlt, & + maero, mpond, & + l_stop, & + istop, jstop, & + aredistn, vredistn) + + if (l_stop) return + + !----------------------------------------------------------------- + ! Make sure the new area = 1. If not (because the closing + ! and opening rates were reduced above), prepare to ridge again + ! with new rates. + !----------------------------------------------------------------- + + call asum_ridging (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, aice0, & + asum) + + call ridge_check (nx_block, ny_block, & + icells, indxi, indxj, & + dt, & + asum, closing_net, & + divu_adv, opning, & + iterate_ridging) + + !----------------------------------------------------------------- + ! If done, exit. If not, prepare to ridge again. + !----------------------------------------------------------------- + +!jd if (iterate_ridging) then +!jd write(nu_diag,*) 'Repeat ridging, niter =', niter +!jd else +!jd exit +!jd endif + if (.not. iterate_ridging) exit + + if (niter == nitermax) then + write(nu_diag,*) ' ' + write(nu_diag,*) 'Exceeded max number of ridging iterations' + write(nu_diag,*) 'max =',nitermax + l_stop = .true. + return + endif + + enddo ! niter + + !----------------------------------------------------------------- + ! Compute final values of conserved quantities. + ! Check for conservation (allowing for snow thrown into ocean). + !----------------------------------------------------------------- + + if (l_conservation_check) then + + eicen(:,:,:) = c0 + esnon(:,:,:) = c0 + sicen(:,:,:) = c0 + vbrin(:,:,:) = c0 + + do n = 1, ncat + do k = 1, nilyr + do j = 1, ny_block + do i = 1, nx_block + eicen(i,j,n) = eicen(i,j,n) + trcrn(i,j,nt_qice+k-1,n) & + * vicen(i,j,n)/real(nilyr,kind=dbl_kind) + enddo + enddo + enddo + do k = 1, nslyr + do j = 1, ny_block + do i = 1, nx_block + esnon(i,j,n) = esnon(i,j,n) + trcrn(i,j,nt_qsno+k-1,n) & + * vsnon(i,j,n)/real(nslyr,kind=dbl_kind) + enddo + enddo + enddo + + + do j = 1, ny_block + do i = 1, nx_block + vbrin(i,j,n) = vicen(i,j,n) + if (tr_brine) vbrin(i,j,n) = trcrn(i,j,nt_fbri,n) * vbrin(i,j,n) + enddo + enddo + + + do k = 1, nilyr + do j = 1, ny_block + do i = 1, nx_block + sicen(i,j,n) = sicen(i,j,n) + trcrn(i,j,nt_sice+k-1,n) & + * vicen(i,j,n)/real(nilyr,kind=dbl_kind) + enddo + enddo + enddo + + enddo + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vicen, vice_final) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vsnon, vsno_final) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + eicen, eice_final) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + esnon, esno_final) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + sicen, Shi_final) + + call column_sum (nx_block, ny_block, & + icells, indxi, indxj, & + ncat, & + vbrin, vbri_final) + + do ij = 1, icells + vsno_final(ij) = vsno_final(ij) + msnow_mlt(ij)/rhos + esno_final(ij) = esno_final(ij) + esnow_mlt(ij) + enddo + + fieldid = 'vice, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + vice_init, vice_final, & + puny, l_stop, & + istop, jstop) + if (l_stop) return + + fieldid = 'vsno, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + vsno_init, vsno_final, & + puny, l_stop, & + istop, jstop) + if (l_stop) return + + fieldid = 'eice, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + eice_init, eice_final, & + puny*Lfresh*rhoi, & + l_stop, & + istop, jstop) + if (l_stop) return + + fieldid = 'esno, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + esno_init, esno_final, & + puny*Lfresh*rhos, & + l_stop, & + istop, jstop) + if (l_stop) return + + fieldid = 'sice, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + Shi_init, Shi_final, & + puny, l_stop, & + istop, jstop) + if (l_stop) return + + fieldid = 'vbrin, ridging' + call column_conservation_check (nx_block, ny_block, & + icells, indxi, indxj, & + fieldid, & + vbri_init, vbri_final, & + puny*c10, l_stop, & + istop, jstop) + if (l_stop) return + endif ! l_conservation_check + + !----------------------------------------------------------------- + ! Compute ridging diagnostics. + !----------------------------------------------------------------- + + dti = c1/dt + + if (present(dardg1dt)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dardg1dt(i,j) = ardg1(ij)*dti + enddo + endif + if (present(dardg2dt)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dardg2dt(i,j) = ardg2(ij)*dti + enddo + endif + if (present(dvirdgdt)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dvirdgdt(i,j) = virdg(ij)*dti + enddo + endif + if (present(opening)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + opening(i,j) = aopen(ij)*dti + enddo + endif + if (present(dardg1ndt)) then + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dardg1ndt(i,j,n) = ardg1n(ij,n)*dti + enddo + enddo + endif + if (present(dardg2ndt)) then + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dardg2ndt(i,j,n) = ardg2n(ij,n)*dti + enddo + enddo + endif + if (present(dvirdgndt)) then + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + dvirdgndt(i,j,n) = virdgn(ij,n)*dti + enddo + enddo + endif + if (present(araftn)) then + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + araftn(i,j,n) = mraftn(ij,n)*ardg2n(ij,n) +! araftn(i,j,n) = mraftn(ij,n)*ardg1n(ij,n)*p5 + enddo + enddo + endif + if (present(vraftn)) then + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + vraftn(i,j,n) = mraftn(ij,n)*virdgn(ij,n) + enddo + enddo + endif + + !----------------------------------------------------------------- + ! Update fresh water and heat fluxes due to snow melt. + !----------------------------------------------------------------- + + ! use thermodynamic time step (ndtd*dt here) to average properly + dti = c1/(ndtd*dt) + + if (present(fresh)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + fresh(i,j) = fresh(i,j) + msnow_mlt(ij)*dti + enddo + endif + if (present(fhocn)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + fhocn(i,j) = fhocn(i,j) + esnow_mlt(ij)*dti + enddo + endif + if (present(faero_ocn)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + faero_ocn(i,j,:) = faero_ocn(i,j,:) + maero(ij,:)*dti + enddo + endif + if (present(fpond)) then + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + fpond(i,j) = fpond(i,j) - mpond(ij) ! units change later + enddo + endif + + !----------------------------------------------------------------- + ! Check for fractional ice area > 1. + !----------------------------------------------------------------- + + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (abs(asum(ij) - c1) > puny) then + l_stop = .true. + istop = i + jstop = j + + write(nu_diag,*) ' ' + write(nu_diag,*) 'Ridging error: total area > 1' + write(nu_diag,*) 'i, j, area:', i, j, asum(ij) + write(nu_diag,*) 'n, aicen:' + write(nu_diag,*) 0, aice0(i,j) + do n = 1, ncat + write(nu_diag,*) n, aicen(i,j,n) + enddo + return + endif + enddo + + end subroutine ridge_ice + +!======================================================================= + +! Find the total area of ice plus open water in each grid cell. +! +! This is similar to the aggregate_area subroutine except that the +! total area can be greater than 1, so the open water area is +! included in the sum instead of being computed as a residual. +! +! author: William H. Lipscomb, LANL + + subroutine asum_ridging (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, aice0, & + asum) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells ! number of cells with ice present + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(in) :: & + aicen ! concentration of ice in each category + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + aice0 ! concentration of open water + + real (kind=dbl_kind), dimension (icells), intent(out):: & + asum ! sum of ice and open water area + + ! local variables + + integer (kind=int_kind) :: & + i, j, n, & + ij ! horizontal index, combines i and j loops + + !----------------------------------------------------------------- + ! open water + !----------------------------------------------------------------- + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + asum(ij) = aice0(i,j) + enddo + + !----------------------------------------------------------------- + ! ice categories + !----------------------------------------------------------------- + + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + asum(ij) = asum(ij) + aicen(i,j,n) + enddo + enddo + + end subroutine asum_ridging + +!======================================================================= + +! Initialize arrays, compute area of closing and opening +! +! author: William H. Lipscomb, LANL + + subroutine ridge_prep (nx_block, ny_block, & + icells, indxi, indxj, & + dt, & + rdg_conv, rdg_shear, & + asum, closing_net, & + divu_adv, opning) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells ! number of cells with ice present + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + + real (kind=dbl_kind), intent(in) :: & + dt ! time step (s) + + real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & + rdg_conv, & ! normalized energy dissipation due to convergence (1/s) + rdg_shear ! normalized energy dissipation due to shear (1/s) + + real (kind=dbl_kind), dimension(icells), & + intent(inout):: & + asum ! sum of ice and open water area + + real (kind=dbl_kind), dimension(icells), & + intent(out):: & + closing_net, & ! net rate at which area is removed (1/s) + divu_adv , & ! divu as implied by transport scheme (1/s) + opning ! rate of opening due to divergence/shear + + ! local variables + + real (kind=dbl_kind), parameter :: & + big = 1.0e+8_dbl_kind + + integer (kind=int_kind) :: & + i,j, & ! horizontal indices + ij ! horizontal index, combines i and j loops + + ! Set hin_max(ncat) to a big value to ensure that all ridged ice + ! is thinner than hin_max(ncat). + hin_max(ncat) = big + + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + !----------------------------------------------------------------- + ! Compute the net rate of closing due to convergence + ! and shear, based on Flato and Hibler (1995). + ! + ! For the elliptical yield curve: + ! rdg_conv = -min (divu, 0) + ! rdg_shear = (1/2) * (Delta - abs(divu)) + ! Note that the shear term also accounts for divergence. + ! + ! The energy dissipation rate is equal to the net closing rate + ! times the ice strength. + ! + ! NOTE: The NET closing rate is equal to the rate that open water + ! area is removed, plus the rate at which ice area is removed by + ! ridging, minus the rate at which area is added in new ridges. + ! The GROSS closing rate is equal to the first two terms (open + ! water closing and thin ice ridging) without the third term + ! (thick, newly ridged ice). + ! + ! rdg_conv is calculated differently in EAP (update_ice_rdg) and + ! represents closing_net directly. In that case, rdg_shear=0. + !----------------------------------------------------------------- + + closing_net(ij) = Cs*rdg_shear(i,j) + rdg_conv(i,j) + + !----------------------------------------------------------------- + ! Compute divu_adv, the divergence rate given by the transport/ + ! advection scheme, which may not be equal to divu as computed + ! from the velocity field. + ! + ! If divu_adv < 0, make sure the closing rate is large enough + ! to give asum = 1.0 after ridging. + !----------------------------------------------------------------- + + divu_adv(ij) = (c1-asum(ij)) / dt + + if (divu_adv(ij) < c0) & + closing_net(ij) = max(closing_net(ij), -divu_adv(ij)) + + !----------------------------------------------------------------- + ! Compute the (non-negative) opening rate that will give + ! asum = 1.0 after ridging. + !----------------------------------------------------------------- + opning(ij) = closing_net(ij) + divu_adv(ij) + + enddo + + end subroutine ridge_prep + +!======================================================================= + +! Compute the thickness distribution of the ice and open water +! participating in ridging and of the resulting ridges. +! +! This version includes new options for ridging participation and +! redistribution. +! The new participation scheme (krdg_partic = 1) improves stability +! by increasing the time scale for large changes in ice strength. +! The new exponential redistribution function (krdg_redist = 1) improves +! agreement between ITDs of modeled and observed ridges. +! +! author: William H. Lipscomb, LANL +! +! 2006: Changed subroutine name to ridge_itd +! Added new options for ridging participation and redistribution. + + subroutine ridge_itd (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, vicen, & + aice0, & + aksum, apartic, & + hrmin, hrmax, & + hrexp, krdg, & + aparticn, krdgn, mraft) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells ! number of cells with ice present + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(in) :: & + aicen , & ! concentration of ice + vicen ! volume per unit area of ice (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + aice0 ! concentration of open water + + real (kind=dbl_kind), dimension (icells), intent(out):: & + aksum ! ratio of area removed to area ridged + + real (kind=dbl_kind), dimension (icells,0:ncat), & + intent(out) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (icells,ncat), & + intent(out) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg ! mean ridge thickness/thickness of ridging ice + + ! diagnostic, category values + real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), & + intent(out), optional :: & + aparticn , & ! participation function + krdgn ! mean ridge thickness/thickness of ridging ice + + real (kind=dbl_kind), dimension (icells,ncat), & + intent(out), optional :: & + mraft ! rafting ice mask + + ! local variables + + integer (kind=int_kind) :: & + i,j , & ! horizontal indices + n , & ! thickness category index + ij ! horizontal index, combines i and j loops + + real (kind=dbl_kind), parameter :: & + Gstari = c1/Gstar, & + astari = c1/astar + + real (kind=dbl_kind), dimension(icells,-1:ncat) :: & + Gsum ! Gsum(n) = sum of areas in categories 0 to n + + real (kind=dbl_kind), dimension(icells) :: & + work ! temporary work array + + real (kind=dbl_kind) :: & + hi , & ! ice thickness for each cat (m) + hrmean , & ! mean ridge thickness (m) + xtmp ! temporary variable + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + + do ij = 1, icells + Gsum (ij,-1) = c0 ! by definition + Gsum (ij,0) = c1 ! to avoid divzero below + apartic(ij,0) = c0 + enddo + + do n = 1, ncat + do ij = 1, icells + Gsum (ij,n) = c1 ! to avoid divzero below + apartic(ij,n) = c0 + hrmin (ij,n) = c0 + hrmax (ij,n) = c0 + hrexp (ij,n) = c0 + krdg (ij,n) = c1 + enddo + enddo + + !----------------------------------------------------------------- + ! Compute the thickness distribution of ice participating in ridging. + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! First compute the cumulative thickness distribution function Gsum, + ! where Gsum(n) is the fractional area in categories 0 to n. + ! Ignore categories with very small areas. + !----------------------------------------------------------------- + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aice0(i,j) > puny) then + Gsum(ij,0) = aice0(i,j) + else + Gsum(ij,0) = Gsum(ij,-1) + endif + enddo + + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aicen(i,j,n) > puny) then + Gsum(ij,n) = Gsum(ij,n-1) + aicen(i,j,n) + else + Gsum(ij,n) = Gsum(ij,n-1) + endif + enddo + enddo + + ! normalize + + do ij = 1, icells + work(ij) = c1 / Gsum(ij,ncat) + enddo + do n = 0, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + Gsum(ij,n) = Gsum(ij,n) * work(ij) + enddo + enddo + + !----------------------------------------------------------------- + ! Compute the participation function apartic; this is analogous to + ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). + ! + ! area lost from category n due to ridging/closing + ! apartic(n) = -------------------------------------------------- + ! total area lost due to ridging/closing + ! + !----------------------------------------------------------------- + + if (krdg_partic == 0) then ! Thornike et al. 1975 formulation + + !----------------------------------------------------------------- + ! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). + ! The expressions for apartic are found by integrating b(h)g(h) between + ! the category boundaries. + !----------------------------------------------------------------- + + do n = 0, ncat + do ij = 1, icells + if (Gsum(ij,n) < Gstar) then + apartic(ij,n) = Gstari*(Gsum(ij,n)-Gsum(ij,n-1)) * & + (c2 - (Gsum(ij,n-1)+Gsum(ij,n))*Gstari) + elseif (Gsum(ij,n-1) < Gstar) then + apartic(ij,n) = Gstari * (Gstar-Gsum(ij,n-1)) * & + (c2 - (Gsum(ij,n-1)+Gstar)*Gstari) + endif + enddo ! ij + enddo ! n + + elseif (krdg_partic==1) then ! exponential dependence on G(h) + + !----------------------------------------------------------------- + ! b(h) = exp(-G(h)/astar) + ! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)]. + ! The expression for apartic is found by integrating b(h)g(h) + ! between the category boundaries. + !----------------------------------------------------------------- + + ! precompute exponential terms using Gsum as work array + + xtmp = c1 / (c1 - exp(-astari)) + + do n = -1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + Gsum(ij,n) = exp(-Gsum(ij,n)*astari) * xtmp + enddo ! ij + enddo ! n + + do n = 0, ncat + do ij = 1, icells + apartic(ij,n) = Gsum(ij,n-1) - Gsum(ij,n) + enddo ! ij + enddo ! n + + endif ! krdg_partic + + !----------------------------------------------------------------- + ! Compute variables related to ITD of ridged ice: + ! + ! krdg = mean ridge thickness / thickness of ridging ice + ! hrmin = min ridge thickness + ! hrmax = max ridge thickness (krdg_redist = 0) + ! hrexp = ridge e-folding scale (krdg_redist = 1) + !---------------------------------------------------------------- + + if (krdg_redist == 0) then ! Hibler 1980 formulation + + !----------------------------------------------------------------- + ! Assume ridged ice is uniformly distributed between hrmin and hrmax. + ! + ! This parameterization is a modified version of Hibler (1980). + ! In the original paper the min ridging thickness is hrmin = 2*hi, + ! and the max thickness is hrmax = 2*sqrt(hi*Hstar). + ! + ! Here the min thickness is hrmin = min(2*hi, hi+maxraft), + ! so thick ridging ice is not required to raft. + ! + !----------------------------------------------------------------- + + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + if (aicen(i,j,n) > puny) then + hi = vicen(i,j,n) / aicen(i,j,n) + hrmin(ij,n) = min(c2*hi, hi + maxraft) + hrmax(ij,n) = c2*sqrt(Hstar*hi) + hrmax(ij,n) = max(hrmax(ij,n), hrmin(ij,n)+puny) + hrmean = p5 * (hrmin(ij,n) + hrmax(ij,n)) + krdg(ij,n) = hrmean / hi + + ! diagnostic rafting mask not implemented + endif + + enddo ! ij + enddo ! n + + else ! krdg_redist = 1; exponential redistribution + + !----------------------------------------------------------------- + ! The ridge ITD is a negative exponential: + ! + ! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin + ! + ! where hrmin is the minimum thickness of ridging ice and + ! hrexp is the e-folding thickness. + ! + ! Here, assume as above that hrmin = min(2*hi, hi+maxraft). + ! That is, the minimum ridge thickness results from rafting, + ! unless the ice is thicker than maxraft. + ! + ! Also, assume that hrexp = mu_rdg*sqrt(hi). + ! The parameter mu_rdg is tuned to give e-folding scales mostly + ! in the range 2-4 m as observed by upward-looking sonar. + ! + ! Values of mu_rdg in the right column give ice strengths + ! roughly equal to values of Hstar in the left column + ! (within ~10 kN/m for typical ITDs): + ! + ! Hstar mu_rdg + ! + ! 25 3.0 + ! 50 4.0 + ! 75 5.0 + ! 100 6.0 + !----------------------------------------------------------------- + + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aicen(i,j,n) > puny) then + hi = vicen(i,j,n) / aicen(i,j,n) + hi = max(hi,puny) + hrmin(ij,n) = min(c2*hi, hi + maxraft) + hrexp(ij,n) = mu_rdg * sqrt(hi) + krdg(ij,n) = (hrmin(ij,n) + hrexp(ij,n)) / hi + + !echmod: check computational efficiency + ! diagnostic rafting mask + if (present(mraft)) then + mraft(ij,n) = max(c0, sign(c1, hi+maxraft-hrmin(ij,n))) + xtmp = mraft(ij,n)*((c2*hi+hrexp(ij,n))/hi - krdg(ij,n)) + mraft(ij,n) = max(c0, sign(c1, puny-abs(xtmp))) + endif + endif + enddo + enddo + + endif ! krdg_redist + + !---------------------------------------------------------------- + ! Compute aksum = net ice area removed / total area participating. + ! For instance, if a unit area of ice with h = 1 participates in + ! ridging to form a ridge with a = 1/3 and h = 3, then + ! aksum = 1 - 1/3 = 2/3. + !---------------------------------------------------------------- + + do ij = 1, icells + aksum(ij) = apartic(ij,0) ! area participating = area removed + enddo + + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + ! area participating > area removed + aksum(ij) = aksum(ij) & + + apartic(ij,n) * (c1 - c1/krdg(ij,n)) + enddo + enddo + + ! diagnostics + if (present(aparticn)) then + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + aparticn(i,j,n) = apartic(ij,n) + enddo + enddo + endif + if (present(krdgn)) then + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + krdgn(i,j,n) = krdg(ij,n) + enddo + enddo + endif + + end subroutine ridge_itd + +!======================================================================= + +! Remove area, volume, and energy from each ridging category +! and add to thicker ice categories. +! +! Tracers: Ridging conserves ice volume and therefore conserves volume +! tracers. It does not conserve ice area, and therefore a portion of area +! tracers are lost (corresponding to the net closing). Area tracers on +! ice that participates in ridging are carried onto the resulting ridged +! ice (except the portion that are lost due to closing). Therefore, +! tracers must be decremented if they are lost to the ocean during ridging +! (e.g. snow, ponds) or if they are being carried only on the level ice +! area. +! +! author: William H. Lipscomb, LANL + + subroutine ridge_shift (nx_block, ny_block, & + icells, indxi, indxj, & + ntrcr, dt, & + aicen, trcrn, & + vicen, vsnon, & + aice0, trcr_depend, & + aksum, apartic, & + hrmin, hrmax, & + hrexp, krdg, & + closing_net, opning, & + ardg1, ardg2, & + virdg, aopen, & + ardg1nn, ardg2nn, & + virdgnn, & + msnow_mlt, esnow_mlt, & + maero, mpond, & + l_stop, & + istop, jstop, & + aredistn, vredistn) + + use ice_state, only: nt_qsno, & + nt_alvl, nt_vlvl, nt_aero, tr_lvl, tr_aero, & + nt_apnd, nt_hpnd, nt_ipnd, tr_pond, & + tr_pond_cesm, tr_pond_lvl, tr_pond_topo, & + nt_fbri + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells , & ! number of cells with ice present + ntrcr ! number of tracers in use + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + real (kind=dbl_kind), intent(in) :: & + dt ! time step (s) + + integer (kind=int_kind), dimension (ntrcr), intent(in) :: & + trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon + + real (kind=dbl_kind), dimension (nx_block,ny_block), & + intent(inout) :: & + aice0 ! concentration of open water + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(inout) :: & + aicen , & ! concentration of ice + vicen , & ! volume per unit area of ice (m) + vsnon ! volume per unit area of snow (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), & + intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), dimension (icells), intent(in) :: & + aksum ! ratio of area removed to area ridged + + real (kind=dbl_kind), dimension (icells,0:ncat), intent(in) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (icells,ncat), intent(in) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg ! mean ridge thickness/thickness of ridging ice + + real (kind=dbl_kind), dimension(icells), intent(inout) :: & + closing_net, & ! net rate at which area is removed (1/s) + opning , & ! rate of opening due to divergence/shear (1/s) + ardg1 , & ! fractional area loss by ridging ice + ardg2 , & ! fractional area gain by new ridges + virdg , & ! ice volume ridged (m) + aopen ! area opened due to divergence/shear + + real (kind=dbl_kind), dimension(icells,ncat), intent(inout) :: & + ardg1nn , & ! area of ice ridged + ardg2nn , & ! area of new ridges + virdgnn ! ridging ice volume + + real (kind=dbl_kind), dimension(icells), intent(inout) :: & + msnow_mlt, & ! mass of snow added to ocean (kg m-2) + esnow_mlt, & ! energy needed to melt snow in ocean (J m-2) + mpond ! mass of pond added to ocean (kg m-2) + + real (kind=dbl_kind), dimension(icells,max_aero), intent(inout) :: & + maero ! aerosol mass added to ocean (kg m-2) + + logical (kind=log_kind), intent(inout) :: & + l_stop ! if true, abort on return + + integer (kind=int_kind), intent(inout) :: & + istop, jstop ! indices of grid cell where model aborts + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(inout), optional :: & + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume + + ! local variables + + integer (kind=int_kind) :: & + i,j , & ! horizontal indices + n, nr , & ! thickness category indices + k , & ! ice layer index + it , & ! tracer index + ij, m , & ! horizontal indices, combine i and j loops + iridge ! number of cells with nonzero ridging + + integer (kind=int_kind), dimension (icells) :: & + indxii, indxjj, & ! compressed indices + indxij ! compressed indices + + real (kind=dbl_kind), dimension (icells,ncat) :: & + aicen_init , & ! ice area before ridging + vicen_init , & ! ice volume before ridging + vsnon_init ! snow volume before ridging + + real (kind=dbl_kind), dimension(icells,ntrcr,ncat) :: & + atrcrn ! aicen*trcrn + + real (kind=dbl_kind), dimension (icells) :: & + closing_gross ! rate at which area removed, not counting + ! area of new ridges + +! ECH note: the following arrays only need be defined on iridge cells + real (kind=dbl_kind), dimension (icells) :: & + afrac , & ! fraction of category area ridged + ardg1n , & ! area of ice ridged + ardg2n , & ! area of new ridges + virdgn , & ! ridging ice volume + vsrdgn , & ! ridging snow volume + dhr , & ! hrmax - hrmin + dhr2 , & ! hrmax^2 - hrmin^2 + farea , & ! fraction of new ridge area going to nr + fvol ! fraction of new ridge volume going to nr + + real (kind=dbl_kind) :: & + esrdgn ! ridging snow energy + + real (kind=dbl_kind) :: & + hi1 , & ! thickness of ridging ice + hexp , & ! ridge e-folding thickness + hL, hR , & ! left and right limits of integration + expL, expR , & ! exponentials involving hL, hR + tmpfac , & ! factor by which opening/closing rates are cut + wk1 ! work variable + + !----------------------------------------------------------------- + ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn + !----------------------------------------------------------------- + + do n = 1, ncat + do it = 1, ntrcr + if (trcr_depend(it) == 0) then ! ice area tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = aicen(i,j,n)*trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 1) then ! ice volume tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = vicen(i,j,n)*trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2) then ! snow volume tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = vsnon(i,j,n)*trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = vicen(i,j,n) & + * trcrn(i,j,nt_fbri,n) & + * trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2+nt_alvl) then ! level ice tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = aicen(i,j,n) & + * trcrn(i,j,nt_alvl,n) & + * trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2+nt_apnd .and. & + (tr_pond_cesm .or. tr_pond_topo)) then ! CESM or topo pond area tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = aicen(i,j,n) & + * trcrn(i,j,nt_apnd,n) & + * trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2+nt_apnd .and. & + tr_pond_lvl) then ! level-ice pond area tracer + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + atrcrn(ij,it,n) = aicen(i,j,n) & + * trcrn(i,j,nt_alvl,n) & + * trcrn(i,j,nt_apnd,n) & + * trcrn(i,j,it,n) + enddo + endif + enddo + enddo + + !----------------------------------------------------------------- + ! Based on the ITD of ridging and ridged ice, convert the net + ! closing rate to a gross closing rate. + ! NOTE: 0 < aksum <= 1 + !----------------------------------------------------------------- + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + closing_gross(ij) = closing_net(ij) / aksum(ij) + + !----------------------------------------------------------------- + ! Reduce the closing rate if more than 100% of the open water + ! would be removed. Reduce the opening rate proportionately. + !----------------------------------------------------------------- + + if (apartic(ij,0) > c0) then + wk1 = apartic(ij,0) * closing_gross(ij) * dt + if (wk1 > aice0(i,j)) then + tmpfac = aice0(i,j) / wk1 + closing_gross(ij) = closing_gross(ij) * tmpfac + opning(ij) = opning(ij) * tmpfac + endif + endif + + enddo ! ij + + !----------------------------------------------------------------- + ! Reduce the closing rate if more than 100% of any ice category + ! would be removed. Reduce the opening rate proportionately. + !----------------------------------------------------------------- + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0) then + wk1 = apartic(ij,n) * closing_gross(ij) * dt + if (wk1 > aicen(i,j,n)) then + tmpfac = aicen(i,j,n) / wk1 + closing_gross(ij) = closing_gross(ij) * tmpfac + opning(ij) = opning(ij) * tmpfac + endif + endif + + enddo ! ij + enddo ! n + + !----------------------------------------------------------------- + ! Compute change in open water area due to closing and opening. + !----------------------------------------------------------------- + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + aice0(i,j) = aice0(i,j) & + - apartic(ij,0)*closing_gross(ij)*dt & + + opning(ij)*dt + if (aice0(i,j) < -puny) then + l_stop = .true. + istop = i + jstop = j + + write (nu_diag,*) ' ' + write (nu_diag,*) 'Ridging error: aice0 < 0' + write (nu_diag,*) 'i, j, aice0:', i, j, aice0(i,j) + return + + elseif (aice0(i,j) < c0) then ! roundoff error + aice0(i,j) = c0 + endif + + aopen(ij) = opning(ij)*dt ! optional diagnostic + + enddo + + !----------------------------------------------------------------- + ! Save initial state variables + !----------------------------------------------------------------- + + do n = 1, ncat + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + aicen_init(ij,n) = aicen(i,j,n) + vicen_init(ij,n) = vicen(i,j,n) + vsnon_init(ij,n) = vsnon(i,j,n) + enddo + enddo + + !----------------------------------------------------------------- + ! Compute the area, volume, and energy of ice ridging in each + ! category, along with the area of the resulting ridge. + !----------------------------------------------------------------- + + do n = 1, ncat + + !----------------------------------------------------------------- + ! Identify grid cells with nonzero ridging + !----------------------------------------------------------------- + + iridge = 0 + + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aicen_init(ij,n) > puny .and. apartic(ij,n) > c0 & + .and. closing_gross(ij) > c0) then + iridge = iridge + 1 + indxii(iridge) = i + indxjj(iridge) = j + indxij(iridge) = ij + endif + enddo ! ij + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + + !----------------------------------------------------------------- + ! Compute area of ridging ice (ardg1n) and of new ridge (ardg2n). + ! Make sure ridging fraction <=1. (Roundoff errors can give + ! ardg1 slightly greater than aicen.) + !----------------------------------------------------------------- + + ardg1n(ij) = apartic(m,n)*closing_gross(m)*dt + + if (ardg1n(ij) > aicen_init(m,n) + puny) then + l_stop = .true. + istop = i + jstop = j + + write (nu_diag,*) ' ' + write (nu_diag,*) 'Ridging error: ardg > aicen' + write (nu_diag,*) 'i, j, n:', i, j, n + write (nu_diag,*) 'ardg, aicen:', & + ardg1n(ij), aicen_init(m,n) + return + else + ardg1n(ij) = min(aicen_init(m,n), ardg1n(ij)) + endif + + ardg2n(ij) = ardg1n(ij) / krdg(m,n) + afrac(ij) = ardg1n(ij) / aicen_init(m,n) + + !----------------------------------------------------------------- + ! Subtract area, volume, and energy from ridging category n. + ! Note: Tracer values are unchanged. + !----------------------------------------------------------------- + + virdgn(ij) = vicen_init(m,n) * afrac(ij) + vsrdgn(ij) = vsnon_init(m,n) * afrac(ij) + + aicen(i,j,n) = aicen(i,j,n) - ardg1n(ij) + vicen(i,j,n) = vicen(i,j,n) - virdgn(ij) + vsnon(i,j,n) = vsnon(i,j,n) - vsrdgn(ij) + + !----------------------------------------------------------------- + ! Increment ridging diagnostics + !----------------------------------------------------------------- + + ardg1(m) = ardg1(m) + ardg1n(ij) + ardg2(m) = ardg2(m) + ardg2n(ij) + virdg(m) = virdg(m) + virdgn(ij) + + ardg1nn(m,n) = ardg1n(ij) + ardg2nn(m,n) = ardg2n(ij) + virdgnn(m,n) = virdgn(ij) + + !----------------------------------------------------------------- + ! Place part of the snow and tracer lost by ridging into the ocean. + !----------------------------------------------------------------- + + msnow_mlt(m) = msnow_mlt(m) + rhos*vsrdgn(ij)*(c1-fsnowrdg) + + if (tr_aero) then + do it = 1, n_aero + maero(m,it) = maero(m,it) & + + vsrdgn(ij)*(c1-fsnowrdg) & + *(trcrn(i,j,nt_aero +4*(it-1),n) & + + trcrn(i,j,nt_aero+1+4*(it-1),n)) + enddo + endif + + if (tr_pond_topo) then + mpond(m) = mpond(m) + ardg1n(ij) & + * trcrn(i,j,nt_apnd,n) & + * trcrn(i,j,nt_hpnd,n) + endif + + !----------------------------------------------------------------- + ! Compute quantities used to apportion ice among categories + ! in the nr loop below + !----------------------------------------------------------------- + + dhr(ij) = hrmax(ij,n) - hrmin(m,n) + dhr2(ij) = hrmax(ij,n) * hrmax(ij,n) - hrmin(m,n) * hrmin(m,n) + + enddo ! ij + + !----------------------------------------------------------------- + ! Increment energy needed to melt snow in ocean. + ! Note that esnow_mlt < 0; the ocean must cool to melt snow. + !----------------------------------------------------------------- + + do k = 1, nslyr +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + esrdgn = vsrdgn(ij) * trcrn(i,j,nt_qsno+k-1,n) & + / real(nslyr,kind=dbl_kind) + esnow_mlt(m) = esnow_mlt(m) + esrdgn*(c1-fsnowrdg) + enddo + enddo + + !----------------------------------------------------------------- + ! Subtract area- and volume-weighted tracers from category n. + !----------------------------------------------------------------- + + do it = 1, ntrcr + if (trcr_depend(it) == 0) then ! ice area tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - ardg1n(ij)*trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 1) then ! ice volume tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - virdgn(ij)*trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 2) then ! snow volume tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - vsrdgn(ij)*trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 2+nt_alvl) then ! level ice tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - ardg1n(ij)*trcrn(i,j,nt_alvl,n)*trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 2+nt_apnd .and. & + (tr_pond_cesm .or. tr_pond_topo)) then ! CESM or topo pond area tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - ardg1n(ij)*trcrn(i,j,nt_apnd,n)*trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 2+nt_apnd .and. & + tr_pond_lvl) then ! level-ice pond area tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - ardg1n(ij) & + * trcrn(i,j,nt_alvl,n) & + * trcrn(i,j,nt_apnd,n) & + * trcrn(i,j,it,n) + enddo + + elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,n) = atrcrn(m,it,n) & + - virdgn(ij)*trcrn(i,j,it,n)*trcrn(i,j,nt_fbri,n) + enddo + endif ! trcr_depend + enddo ! ntrcr + + !----------------------------------------------------------------- + ! Add area, volume, and energy of new ridge to each category nr. + !----------------------------------------------------------------- + + do nr = 1, ncat + + if (krdg_redist == 0) then ! Hibler 1980 formulation + + do ij = 1, iridge + m = indxij(ij) + + !----------------------------------------------------------------- + ! Compute the fraction of ridged ice area and volume going to + ! thickness category nr. + !----------------------------------------------------------------- + + if (hrmin(m,n) >= hin_max(nr) .or. & + hrmax(ij,n) <= hin_max(nr-1)) then + hL = c0 + hR = c0 + else + hL = max (hrmin(m,n), hin_max(nr-1)) + hR = min (hrmax(ij,n), hin_max(nr)) + endif + + farea(ij) = (hR-hL) / dhr(ij) + fvol (ij) = (hR*hR - hL*hL) / dhr2(ij) + + enddo ! ij + + else ! krdg_redist = 1; 2005 exponential formulation + + !----------------------------------------------------------------- + ! Compute the fraction of ridged ice area and volume going to + ! thickness category nr. + !----------------------------------------------------------------- + + if (nr < ncat) then + + do ij = 1, iridge + m = indxij(ij) + + hi1 = hrmin(m,n) + hexp = hrexp(m,n) + + if (hi1 >= hin_max(nr)) then + farea(ij) = c0 + fvol (ij) = c0 + else + hL = max (hi1, hin_max(nr-1)) + hR = hin_max(nr) + expL = exp(-(hL-hi1)/hexp) + expR = exp(-(hR-hi1)/hexp) + farea(ij) = expL - expR + fvol (ij) = ((hL + hexp)*expL & + - (hR + hexp)*expR) / (hi1 + hexp) + endif + enddo ! ij + + else ! nr = ncat + + do ij = 1, iridge + m = indxij(ij) + + hi1 = hrmin(m,n) + hexp = hrexp(m,n) + + hL = max (hi1, hin_max(nr-1)) + expL = exp(-(hL-hi1)/hexp) + farea(ij) = expL + fvol (ij) = (hL + hexp)*expL / (hi1 + hexp) + enddo + + endif ! nr < ncat + + ! diagnostics + if (n ==1) then ! only for thinnest ridging ice + if (present(aredistn)) then + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + aredistn(i,j,nr) = farea(ij)*ardg2n(ij) + enddo + endif + if (present(vredistn)) then + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + vredistn(i,j,nr) = fvol(ij)*virdgn(ij) + enddo + endif + endif + + endif ! krdg_redist + + !----------------------------------------------------------------- + ! Transfer ice area, ice volume, and snow volume to category nr. + !----------------------------------------------------------------- + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + aicen(i,j,nr) = aicen(i,j,nr) + farea(ij)*ardg2n(ij) + vicen(i,j,nr) = vicen(i,j,nr) + fvol(ij) *virdgn(ij) + vsnon(i,j,nr) = vsnon(i,j,nr) & + + fvol(ij)*vsrdgn(ij)*fsnowrdg + enddo + + !----------------------------------------------------------------- + ! Transfer area-weighted and volume-weighted tracers to category nr. + ! Note: The global sum aicen*trcrn of ice area tracers + ! (trcr_depend = 0) is not conserved by ridging. + ! However, ridging conserves the global sum of volume + ! tracers (trcr_depend = 1 or 2). + ! Tracers associated with level ice, or that are otherwise lost + ! from ridging ice, are not transferred. + ! We assume that all pond water is lost from ridging ice. + !----------------------------------------------------------------- + + do it = 1, ntrcr + if (trcr_depend(it) == 0) then ! ice area tracer + if (it /= nt_alvl) then +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,nr) = atrcrn(m,it,nr) & + + farea(ij)*ardg2n(ij)*trcrn(i,j,it,n) + enddo + endif + elseif (trcr_depend(it) == 1) then ! ice volume tracer + if (it /= nt_vlvl) then +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,nr) = atrcrn(m,it,nr) & + + fvol(ij)*virdgn(ij)*trcrn(i,j,it,n) + enddo + endif + elseif (trcr_depend(it) == 2) then ! snow volume tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,nr) = atrcrn(m,it,nr) & + + fvol(ij)*vsrdgn(ij)*fsnowrdg*trcrn(i,j,it,n) + enddo + elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, iridge + i = indxii(ij) + j = indxjj(ij) + m = indxij(ij) + atrcrn(m,it,nr) = atrcrn(m,it,nr) & ! + + fvol(ij)*virdgn(ij)*trcrn(i,j,nt_fbri,n)*trcrn(i,j,it,n) + enddo + endif ! trcr_depend + enddo ! ntrcr + + enddo ! nr (new ridges) + enddo ! n (ridging categories) + + !----------------------------------------------------------------- + ! Compute new tracers + !----------------------------------------------------------------- + + do n = 1, ncat + call compute_tracers (nx_block, ny_block, & + icells, indxi, indxj, & + ntrcr, trcr_depend, & + atrcrn(:,:,n), aicen(:,:, n), & + vicen (:,:, n), vsnon(:,:, n), & + trcrn(:,:,:,n)) + enddo + + end subroutine ridge_shift + +!======================================================================= + +! Make sure ice area <=1. If not, prepare to repeat ridging. +! +! authors William H. Lipscomb, LANL + + subroutine ridge_check (nx_block, ny_block, & + icells, indxi, indxj, & + dt, & + asum, closing_net, & + divu_adv, opning, & + iterate_ridging) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icells ! number of cells with ice present + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with ice + + + real (kind=dbl_kind), intent(in) :: & + dt ! time step (s) + + real (kind=dbl_kind), dimension (icells), intent(in) :: & + asum ! sum of ice and open water area + + real (kind=dbl_kind), dimension (icells), & + intent(inout) :: & + closing_net, & ! net rate at which area is removed (1/s) + divu_adv , & ! divu as implied by transport scheme (1/s) + opning ! rate of opening due to divergence/shear + + logical (kind=log_kind), intent(out) :: & + iterate_ridging ! if true, repeat the ridging + + ! local variables + + integer (kind=int_kind) :: & + ij ! horizontal index, combines i and j loops + + iterate_ridging = .false. + + do ij = 1, icells + if (abs(asum(ij) - c1) < puny) then + closing_net(ij) = c0 + opning (ij) = c0 + else + iterate_ridging = .true. + divu_adv(ij) = (c1 - asum(ij)) / dt + closing_net(ij) = max(c0, -divu_adv(ij)) + opning(ij) = max(c0, divu_adv(ij)) + endif + enddo + + end subroutine ridge_check + +!======================================================================= + +! Compute the strength of the ice pack, defined as the energy (J m-2) +! dissipated per unit area removed from the ice pack under compression, +! and assumed proportional to the change in potential energy caused +! by ridging. +! +! See Rothrock (1975) and Hibler (1980). +! +! For simpler strength parameterization, see this reference: +! Hibler, W. D. III, 1979: A dynamic-thermodynamic sea ice model, +! J. Phys. Oceanog., 9, 817-846. +! +! authors: William H. Lipscomb, LANL +! Elizabeth C. Hunke, LANL + + subroutine ice_strength (nx_block, ny_block, & + ilo, ihi, jlo, jhi, & + icells, & + indxi, indxj, & + aice, vice, & + aice0, aicen, & + vicen, strength) + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + ilo,ihi,jlo,jhi ! beg and end of physical domain + + integer (kind=int_kind), intent(in) :: & + icells ! no. of cells where icetmask = 1 + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxi , & ! compressed index in i-direction + indxj ! compressed index in j-direction + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + aice , & ! concentration of ice + vice , & ! volume per unit area of ice (m) + aice0 ! concentration of open water + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & + intent(in) :: & + aicen , & ! concentration of ice + vicen ! volume per unit area of ice (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & + strength ! ice strength (N/m) + + ! local variables + + real (kind=dbl_kind), dimension (icells) :: & + asum , & ! sum of ice and open water area + aksum ! ratio of area removed to area ridged + + real (kind=dbl_kind), dimension (icells,0:ncat) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (icells,ncat) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg ! mean ridge thickness/thickness of ridging ice + + integer (kind=int_kind) :: & + i,j , & ! horizontal indices + n , & ! thickness category index + ij ! horizontal index, combines i and j loops + + real (kind=dbl_kind) :: & + hi , & ! ice thickness (m) + h2rdg , & ! mean value of h^2 for new ridge + dh2rdg ! change in mean value of h^2 per unit area + ! consumed by ridging + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + + strength(:,:) = c0 + + if (kstrength == 1) then ! Rothrock '75 formulation + + !----------------------------------------------------------------- + ! Compute thickness distribution of ridging and ridged ice. + !----------------------------------------------------------------- + + call asum_ridging (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, aice0, & + asum) + + call ridge_itd (nx_block, ny_block, & + icells, indxi, indxj, & + aicen, vicen, & + aice0, & + aksum, apartic, & + hrmin, hrmax, & + hrexp, krdg) + + !----------------------------------------------------------------- + ! Compute ice strength based on change in potential energy, + ! as in Rothrock (1975) + !----------------------------------------------------------------- + + if (krdg_redist==0) then ! Hibler 1980 formulation + + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0)then + hi = vicen(i,j,n) / aicen(i,j,n) + h2rdg = p333 * (hrmax(ij,n)**3 - hrmin(ij,n)**3) & + / (hrmax(ij,n) - hrmin(ij,n)) + dh2rdg = -hi*hi + h2rdg/krdg(ij,n) + strength(i,j) = strength(i,j) & + + apartic(ij,n) * dh2rdg + endif ! aicen > puny + enddo ! ij + enddo ! n + + elseif (krdg_redist==1) then ! exponential formulation + + do n = 1, ncat +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0)then + hi = vicen(i,j,n) / aicen(i,j,n) + h2rdg = hrmin(ij,n)*hrmin(ij,n) & + + c2*hrmin(ij,n)*hrexp(ij,n) & + + c2*hrexp(ij,n)*hrexp(ij,n) + dh2rdg = -hi*hi + h2rdg/krdg(ij,n) + strength(i,j) = strength(i,j) & + + apartic(ij,n) * dh2rdg + endif + enddo ! ij + enddo ! n + + endif ! krdg_redist + +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + strength(i,j) = Cf * Cp * strength(i,j) / aksum(ij) + ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) + ! Cf accounts for frictional dissipation + enddo ! ij + + else ! kstrength /= 1: Hibler (1979) form + + !----------------------------------------------------------------- + ! Compute ice strength as in Hibler (1979) + !----------------------------------------------------------------- + + do j = jlo, jhi + do i = ilo, ihi + strength(i,j) = Pstar*vice(i,j)*exp(-Cstar*(c1-aice(i,j))) + enddo ! j + enddo ! i + + endif ! kstrength + + end subroutine ice_strength + +!======================================================================= + + end module ice_mechred + +!======================================================================= diff --git a/apps/common/modified_src/roms-3.6/frazil_ice_prod_mod.F b/apps/common/modified_src/roms-3.6/frazil_ice_prod_mod.F index 6d70bad..68193f3 100644 --- a/apps/common/modified_src/roms-3.6/frazil_ice_prod_mod.F +++ b/apps/common/modified_src/roms-3.6/frazil_ice_prod_mod.F @@ -5,20 +5,18 @@ MODULE frazil_ice_prod_mod ! !jd- MET.no ! -! This routine computes the frazil ice growth in the water when the -! water temperature gets below freezing by taking account of the energy -! difference between under-cooled water, and water at the freezing point. -! Temperature is then set back to the freezing-temperature. -! This version takes implicitly account of the change in salinity when -! frazil-ice is formed, and the change in salinity is estimated before the -! Available energy for freezing is calculated, and stored. +! Frazil ice is modelled as the energi deficient required to increase the +! water temperature to the freezing-point. We assume that the ice model +! uses all this energy to produce ice, but do no assumption on the amount, +! final temperature or salinity. Salinity changes are therefore calculated +! in the ice model and delivered back to ocean as a surface flux. +! We allow frazil ice energy from deeper layer to be used in cool upper layers +! if their temperature are above freezing. +! We restrict the calculation of frazil ice to the upper Z_R_MAX meter to +! get clear of most of the problematic spurious watermasses generated by +! advection and steep topography. ! ! No known reference at present. Based on calculation by Jens Debernard -! inspired by Steele et al. (1989). JPO, 19, 139-147, -! but without taking account change in ocean enthaply due to mass -! exchange between ocean model and ice model. The present -! version of ROMS, and the coupling is not consistent with a exchange -! of mass. ! ! !======================================================================= ! @@ -36,7 +34,9 @@ real(r8) function t_freeze(s1,z1) ! t_freeze(s1,z1) = -0.0575*s1 + 1.710523d-3*sqrt(s1)**3 ! & - 2.154996d-4*s1*s1 + 0.000753*z1 ! Freezing temperature (Steele et al. 1989) - t_freeze = -0.0543*s1 + 0.000759*z1 +! t_freeze = -0.0543*s1 + 0.000759*z1 +! Freezing temperature, linear in salt, no depth dependence + t_freeze = -0.054_r8*s1 return end function t_freeze @@ -122,7 +122,7 @@ subroutine frazil_ice_prod_tile (ng, tile, & real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) real(r8), intent(out) :: qfraz(LBi:,LBj:) - real(r8), intent(out) :: qfraz_accum(LBi:,LBj:) + real(r8), intent(inout) :: qfraz_accum(LBi:,LBj:) # else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) @@ -134,7 +134,7 @@ subroutine frazil_ice_prod_tile (ng, tile, & real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) real(r8), intent(out) :: qfraz(LBi:UBi,LBj:UBj) - real(r8), intent(out) :: qfraz_accum(LBi:UBi,LBj:UBj) + real(r8), intent(inout) :: qfraz_accum(LBi:UBi,LBj:UBj) # endif ! ! Local variable definitions @@ -144,12 +144,16 @@ subroutine frazil_ice_prod_tile (ng, tile, & real(r8), parameter :: Lh_ice = 334000.0_r8 + real(r8), parameter :: z_r_max= -50.0_r8 + !jd real(r8), parameter :: Lhat = 79.2_r8 !jd real(r8), parameter :: r = 0.5_r8 real(r8) :: rhocp real(r8) :: cpw_hat - real(r8) :: tsuper ! jd Amount of supercooling + real(r8) :: qfraz_prod ! Energy to frazil in this layer (T < T_freeze) + real(r8) :: meltpot ! Energy available to melt frazil in layer (T>T_freeze) + real(r8) :: meltheat ! Energy used to melt frazil in layer real(r8) :: t_fr, Sold @@ -161,39 +165,44 @@ subroutine frazil_ice_prod_tile (ng, tile, & DO j=Jstr,Jend DO i=Istr,Iend qfraz(i,j) = 0.0_r8 - ENDDO - ENDDO - DO j=Jstr,Jend - DO i=Istr,Iend - DO k=1,N(ng) # ifdef MASKING IF (rmask(i,j) .ne. 0.0_r8) THEN # endif # ifdef WET_DRY IF (rmask_wet(i,j) .ne. 0.0_r8) THEN # endif - -!jd To ensure S >= 0, S<0 is unphysical and makes this routine unstable + DO k=1,N(ng) + if (z_r(i,j,k) < z_r_max) cycle +!jd To ensure S >= 0, S<0 is unphysical Sold = max(0.0_r8,t(i,j,k,nnew,isalt)) t_fr = t_freeze(Sold,z_r(i,j,k)) - tsuper= max(t_fr - t(i,j,k,nnew,itemp),0.0_r8) + qfraz_prod= rhocp*Hz(i,j,k) * & + & max(t_fr - t(i,j,k,nnew,itemp),0.0_r8) + + meltpot= rhocp*Hz(i,j,k) * & + & min(t_fr - t(i,j,k,nnew,itemp),0.0_r8) + + +! Add possible frazil ice production + qfraz(i,j) = qfraz(i,j) + qfraz_prod -!jd t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) * & -!jd & (1.0_r8 + cpw_hat*tsuper) +! Allow warm water to melt ice produced in deeper layers - t(i,j,k,nnew,isalt) = Sold * (1.0_r8 + cpw_hat*tsuper) + meltheat= max (meltpot, -qfraz(i,j) ) ! <= 0 - qfraz(i,j) = qfraz(i,j) + rhocp*tsuper*Hz(i,j,k) + qfraz(i,j) = qfraz(i,j) + meltheat - t(i,j,k,nnew,itemp) = max(t(i,j,k,nnew,itemp), & - & t_freeze(t(i,j,k,nnew,isalt),z_r(i,j,k))) +! Adjust temperature in accordance with change in qfraz + t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + & + & (qfraz_prod + meltheat )/(rhocp*Hz(i,j,k)) + + END DO # ifdef WET_DRY END IF # endif # ifdef MASKING END IF # endif - END DO qfraz(i,j) = qfraz(i,j)/dt(ng) IF (qfraz(i,j) .lt. 0.0_r8) THEN print *, 'trouble in frazil_ice_mod', i, j, & diff --git a/apps/common/modified_src/roms-3.6/mod_coupler.F b/apps/common/modified_src/roms-3.6/mod_coupler.F index a3c3d42..ff6a631 100644 --- a/apps/common/modified_src/roms-3.6/mod_coupler.F +++ b/apps/common/modified_src/roms-3.6/mod_coupler.F @@ -195,7 +195,7 @@ SUBROUTINE allocate_coupler (Nnodes) ! ! Open input coupling variable information file. ! - WRITE (stdout,*) ' Opening coupeling variable information file: ', CPLname +!jd WRITE (stdout,*) ' Opening coupeling variable information file: ', CPLname OPEN (inp, FILE=TRIM(CPLname), FORM='formatted', STATUS='old', & & ERR=10) GO TO 20 @@ -206,7 +206,7 @@ SUBROUTINE allocate_coupler (Nnodes) ! Read in variable information. Ignore blank and comment [char(33)=!] ! input lines. ! - WRITE (stdout,*) ' Reading coupling variable information file' +!jd WRITE (stdout,*) ' Reading coupling variable information file' varid=0 DO WHILE (.TRUE.) READ (inp,*,ERR=30,END=40) Vinfo(1) @@ -278,7 +278,7 @@ SUBROUTINE allocate_coupler (Nnodes) ! ! Allocate IDs structures. ! - WRITE (stdout,*) ' Allocate IDs structures' +!jd WRITE (stdout,*) ' Allocate IDs structures' IF (.not.allocated(ExportID)) THEN allocate ( ExportID(Nmodels) ) DO model=1,Nmodels @@ -343,7 +343,7 @@ SUBROUTINE allocate_coupler (Nnodes) ! ! Allocate structure. ! - WRITE (stdout,*) ' Assign processors to coupled models.' +!jd WRITE (stdout,*) ' Assign processors to coupled models.' IF (.not.allocated(pets)) THEN allocate ( pets(Nmodels) ) DO model=1,Nmodels diff --git a/apps/common/modified_src/roms-3.6/ocean_coupler.F b/apps/common/modified_src/roms-3.6/ocean_coupler.F index fbe5509..d2d1182 100644 --- a/apps/common/modified_src/roms-3.6/ocean_coupler.F +++ b/apps/common/modified_src/roms-3.6/ocean_coupler.F @@ -377,10 +377,10 @@ SUBROUTINE ocn2cice_coupling (ng, tile, ncouple) allocate(uwater(LBi:UBi,LBj:UBj)) allocate(vwater(LBi:UBi,LBj:UBj)) -! uw=0.0_r8 -! vw=0.0_r8 -! uwater=0.0_r8 -! vwater=0.0_r8 + uw=0.0_r8 + vw=0.0_r8 + uwater=0.0_r8 + vwater=0.0_r8 IF (Master) THEN @@ -782,7 +782,8 @@ SUBROUTINE ocn2cice_coupling (ng, tile, ncouple) !jd weight over a typical mixed (?)/ near-ice layer. ICE(ng)%qfraz_accum(i,j) = rho0*Cp * & & min(t_freeze(OCEAN(ng)%t(i,j,N(ng),NOUT,isalt),0.0_r8) & - & - OCEAN(ng)%t(i,j,N(ng),NOUT,itemp), 0.0_r8 ) * 5.0_r8 + & - OCEAN(ng)%t(i,j,N(ng),NOUT,itemp), 0.0_r8 ) & + & * 5.0_r8 / (ncouple*dt(ng)) end if end do end do diff --git a/apps/common/origfiles/roms_keyword.in b/apps/common/origfiles/roms_keyword.in index 5c51683..7225279 100644 --- a/apps/common/origfiles/roms_keyword.in +++ b/apps/common/origfiles/roms_keyword.in @@ -584,8 +584,8 @@ Dout(iTvdif) == F F ! vertical diffusion GRDNAME == GRDFILE ININAME == RUNDIR/ocean_ini.nc - CLMNAME == RUNDIR/ocean_clm.nc - BRYNAME == RUNDIR/ocean_bry.nc + CLMNAME == RUNDIR/ocean_1997_clm.nc + BRYNAME == RUNDIR/ocean_1997_bry.nc ! Input forcing NetCDF file name(s). The USER has the option to enter ! several files names per each nested grid. For example, the USER may diff --git a/apps/common/python/ModelRun.py b/apps/common/python/ModelRun.py index 10b761a..f74edbf 100644 --- a/apps/common/python/ModelRun.py +++ b/apps/common/python/ModelRun.py @@ -96,7 +96,7 @@ def _execute_roms_mpi(self,ncpus,infile,debugoption=Constants.NODEBUG,architectu exit(1) else: - os.environ["MPI_BUFS_PER_PROC"] = str(128) +# os.environ["MPI_BUFS_PER_PROC"] = str(128) os.system("mpiexec_mpt -np "+str(ncpus)+" "+executable+" "+infile) else: os.system("mpirun -np "+str(ncpus)+" "+executable+" "+infile) @@ -243,4 +243,4 @@ def _check_starttime(self): print "There seems to be a problem with matching dates in ROMS and CICE. Will exit..." print "ROMS: "+str(roms_ini[self._params.NRREC].day) print "CICE: "+str(cice_rst_day) - exit(1) \ No newline at end of file + exit(1) diff --git a/apps/common/python/Params.py b/apps/common/python/Params.py index 1c8cd03..ddedf53 100644 --- a/apps/common/python/Params.py +++ b/apps/common/python/Params.py @@ -45,7 +45,7 @@ def __init__(self,app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=0,restart=F self.COUPLINGTIME_I2O=3600.0 #self.ROMSINIFILE=self.RUNPATH+"/"+INIFILE # Find restart-time of CICE: - cice_start_step = (start_date-datetime(start_date.year,01,01)).total_seconds()/self.CICEDELTAT + cice_start_step = (start_date-datetime(start_date.year,01,01,00)).total_seconds()/self.CICEDELTAT if restart == True: f = open(self.CICERUNDIR+'/restart/ice.restart_file', 'r') cice_restartfile = f.readline().strip() @@ -68,26 +68,26 @@ def __init__(self,app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=0,restart=F ['NLEVELS',"35"], #Could read from grd-file? ['GRDTHETAS',"6.0d0"], ['GRDTHETAB',"0.1d0"], - ['GRDTCLINE',"30.0d0"], + ['GRDTCLINE',"100.0d0"],# ['XCPU',str(self.XCPU)], ['YCPU',str(self.YCPU)], ['TSTEPS',str(self.FCLEN/self.DELTAT)], ['DELTAT',str(self.DELTAT)], ['RATIO',"20"], #['RATIO',"30"], ['IRESTART',str(self.NRREC)], - ['RSTSTEP',str(24*3600/int(self.DELTAT))], - ['STASTEP',str(1*3600/int(self.DELTAT))], + ['RSTSTEP',str(30*24*3600/int(self.DELTAT))], + ['STASTEP',str(24*3600/int(self.DELTAT))], ['INFOSTEP',str(1*3600/int(self.DELTAT))], ['HISSTEPP',str(1*3600/int(self.DELTAT))], - ['DEFHISSTEP',str(720*3600/int(self.DELTAT))], #if 0; all output in one his-file - ['AVGSTEPP',str(24*3600/int(self.DELTAT))], - ['STARTAVG',"0"], - ['DEFAVGSTEP',str(720*3600/int(self.DELTAT))], #if 0; all output in one avg-file + ['DEFHISSTEP',str(30*24*3600/int(self.DELTAT))], #if 0; all output in one his-file + ['AVGSTEPP',str(1*24*3600/int(self.DELTAT))], + ['STARTAVG',"1"], + ['DEFAVGSTEP',str(30*24*3600/int(self.DELTAT))], #if 0; all output in one avg-file ['STARTTIME',str((start_date-self.TIMEREF).total_seconds()/86400)], ['TIDEREF',str((start_date-self.TIMEREF).total_seconds()/86400)], ['TIMEREF',self.TIMEREF.strftime("%Y%m%d.00")], ['V_TRANS',"2"], - ['V_STRETCH',"1"], + ['V_STRETCH',"2"], ['OBCFAKTOR',"120.0"], ['NUDGZONEWIDTH',"15"], ['GRDFILE',GlobalParams.COMMONPATH+"/grid/A20_grd_openBering.nc"], @@ -138,7 +138,7 @@ def __init__(self,app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=0,restart=F self.CICEDELTAT=600 self.COUPLINGTIME_I2O=600 # Find restart-time of CICE: - cice_start_step = (start_date-datetime(start_date.year,01,01)).total_seconds()/self.CICEDELTAT + cice_start_step = (start_date-datetime(start_date.year,01,01,00)).total_seconds()/self.CICEDELTAT if restart == True: f = open(self.CICERUNDIR+'/restart/ice.restart_file', 'r') cice_restartfile = f.readline().strip() @@ -159,7 +159,7 @@ def __init__(self,app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=0,restart=F ['NLEVELS',"35"], #Could read from grd-file? ['GRDTHETAS',"6.0d0"], ['GRDTHETAB',"0.1d0"], - ['GRDTCLINE',"100.0d0"], + ['GRDTCLINE',"100.0d0"], ['XCPU',str(self.XCPU)], ['YCPU',str(self.YCPU)], ['TSTEPS',str(self.FCLEN/self.DELTAT)], @@ -180,7 +180,7 @@ def __init__(self,app,xcpu,ycpu,start_date,end_date,nrrec=-1,cicecpu=0,restart=F ['TIDEREF',str((start_date-self.TIMEREF).total_seconds()/86400)], ['TIMEREF',self.TIMEREF.strftime("%Y%m%d.00")], ['V_TRANS',"2"], - ['V_STRETCH',"1"], + ['V_STRETCH',"2"], ['OBCFAKTOR',"1.0"], ['NUDGZONEWIDTH',"15"], ['GRDFILE',GlobalParams.COMMONPATH+"/grid/arctic4km_grd.nc"], diff --git a/apps/modules.sh b/apps/modules.sh index 75af4c6..c164851 100755 --- a/apps/modules.sh +++ b/apps/modules.sh @@ -1,13 +1,5 @@ #!/bin/bash source /etc/profile.d/modules.sh -#module load intelcomp/14.0.1 -#module load mpt/2.06 -#module load netcdf/4.1.3 -module load intelcomp mpt netcdf python/2.7.2 -#module unload intelcomp/12.0.5.220 -#module purge -#module load mpt/2.10 intelcomp/15.0.1 netcdf/4.3.2 hdf5/1.8.14-mpi -#module load intelcomp/15.0.1 -#module load netcdf/4.3.2 -#module load hdf5/1.8.14-mpi +module purge +module load mpt/2.06 intelcomp/12.0.5.220 netcdf/4.1.3 python/2.7.2