diff --git a/README.md b/README.md index 457714b8a..45e4a6a1a 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ -[![Travis-CI](https://travis-ci.org/CICE-Consortium/CICE.svg?branch=master)](https://travis-ci.org/CICE-Consortium/CICE) + [![GHActions](https://github.com/CICE-Consortium/CICE/workflows/GHActions/badge.svg)](https://github.com/CICE-Consortium/CICE/actions) -[![Documentation Status](https://readthedocs.org/projects/cice-consortium-cice/badge/?version=master)](http://cice-consortium-cice.readthedocs.io/en/master/?badge=master) +[![Documentation Status](https://readthedocs.org/projects/cice-consortium-cice/badge/?version=main)](http://cice-consortium-cice.readthedocs.io/en/main/?badge=main) [![lcov](https://img.shields.io/endpoint?url=https://apcraig.github.io/coverage.json)](https://apcraig.github.io) diff --git a/cice.setup b/cice.setup index 3cae17d74..60c56e5c2 100755 --- a/cice.setup +++ b/cice.setup @@ -1,5 +1,7 @@ #!/bin/csh -f +#set pd0 = `date -u "+%s%N"` + set ICE_SANDBOX = `pwd` set ICE_VERSION = unknown if (-e cicecore/version.txt) then @@ -480,6 +482,7 @@ else exit -1 endif cp -f ${ICE_SCRIPTS}/tests/report_results.csh ${tsdir} + cp -f ${ICE_SCRIPTS}/tests/create_fails.csh ${tsdir} cp -f ${ICE_SCRIPTS}/tests/poll_queue.csh ${tsdir} cat >! ${tsdir}/suite.submit << EOF0 @@ -823,8 +826,8 @@ EOF # set default test output as failure if (${docase} == 0) then echo "#---" >! test_output - echo "FAIL ${testname_noid} build" >> test_output - echo "FAIL ${testname_noid} run" >> test_output + echo "PEND ${testname_noid} build" >> test_output + echo "PEND ${testname_noid} run" >> test_output endif # from basic script dir to case @@ -919,7 +922,7 @@ EOF echo "ICE_GRID = ${grid} (${ICE_DECOMP_NXGLOB}x${ICE_DECOMP_NYGLOB}) blocksize=${ICE_DECOMP_BLCKX}x${ICE_DECOMP_BLCKY}x${ICE_DECOMP_MXBLCKS}" echo "ICE_DECOMP = ${ICE_DECOMP_DECOMP} ${ICE_DECOMP_DSHAPE}" if ($fbfbcomp != ${spval}) then - echo "ICE_BFBCOMP = ${fbfbcomp}" + echo "ICE_BFBCOMP = ${fbfbcomp}" endif #------------------------------------------------------------ @@ -933,9 +936,21 @@ EOF if (-e ${fimods}) rm ${fimods} if (-e ${fsmods}) rm ${fsmods} + # Use an existing ice_in file from the suite if it exists + # to reduce time spent in parse_namelist + set skip_parse_namelist = spval + if (${dosuite} == 1) then + set iceinfn = ../ice_in_save_${grid}${soptions} + if (-e ${iceinfn}) then + echo "use ${iceinfn}" + cp ${iceinfn} ice_in + set skip_parse_namelist = true + endif + endif + + # Set decomp info in namelist cat >! ${fimods} << EOF1 # cice.setup settings - nprocs = ${task} nx_global = ${ICE_DECOMP_NXGLOB} ny_global = ${ICE_DECOMP_NYGLOB} @@ -964,7 +979,6 @@ EOF1 cat >! ${fsmods} << EOF1 # cice.setup settings - setenv ICE_SANDBOX ${ICE_SANDBOX} setenv ICE_SCRIPTS ${ICE_SCRIPTS} setenv ICE_CASENAME ${casename} @@ -1033,42 +1047,59 @@ EOF1 foreach name (${grid} $setsx) set found = 0 + if (-e ${ICE_SCRIPTS}/options/set_nml.${name}) then cat >> ${fimods} << EOF2 # set_nml.${name} - EOF2 - cat ${ICE_SCRIPTS}/options/set_nml.${name} >> ${fimods} - cat >> ${fimods} << EOF2 - + if ("${skip_parse_namelist}" == "true") then + # need to make sure the decomp info from the set_nml is picked up. each case + # has a slightly different decomp that is independent of the ice_in_save file. + # compute that then overwrite by set_nml as needed. + grep -i "distribution_type" ${ICE_SCRIPTS}/options/set_nml.${name} >> ${fimods} + grep -i "processor_shape" ${ICE_SCRIPTS}/options/set_nml.${name} >> ${fimods} + cat >> ${fimods} << EOF2 +# using saved ice_in EOF2 + else + cat ${ICE_SCRIPTS}/options/set_nml.${name} >> ${fimods} + cat >> ${fimods} << EOF2 +EOF2 + endif echo "adding namelist mods set_nml.${name}" echo "`date` ${0} adding namelist modes set_nml.${name}" >> ${casedir}/README.case set found = 1 endif + if (-e ${ICE_SCRIPTS}/options/set_env.${name}) then cat >> ${fsmods} << EOF2 # set_env.${name} - EOF2 cat ${ICE_SCRIPTS}/options/set_env.${name} >> ${fsmods} cat >> ${fsmods} << EOF2 - EOF2 echo "adding env mods set_env.${name}" echo "`date` ${0} adding namelist modes set_env.${name}" >> ${casedir}/README.case set found = 1 endif + if (${found} == 0) then echo "${0}: ERROR, ${ICE_SCRIPTS}/options/set_[nml,env].${name} not found" exit -1 endif end +#set pd1 = `date -u "+%s%N"` +#@ pdd = ( $pd1 - $pd0 ) / 1000000 +#echo "tcxp b4 parse $pdd" ${casescr}/parse_settings.sh cice.settings ${fsmods} + if ($status != 0) then + echo "${0}: ERROR, parse_settings.sh aborted" + exit -1 + endif ${casescr}/parse_namelist.sh ice_in ${fimods} if ($status != 0) then echo "${0}: ERROR, parse_namelist.sh aborted" @@ -1077,6 +1108,20 @@ EOF2 source ./cice.settings source ./env.${machcomp} -nomodules || exit 2 ${casescr}/parse_namelist_from_env.sh ice_in + if ($status != 0) then + echo "${0}: ERROR, parse_namelist_from_env.sh aborted" + exit -1 + endif +#set pd1 = `date -u "+%s%N"` +#@ pdd = ( $pd1 - $pd0 ) / 1000000 +#echo "tcxp after parse $pdd" + + # Save ice_in in the suite to reduce time spent in parse_namelist + if (${dosuite} == 1) then + if !(-e ${iceinfn}) then + cp ice_in ${iceinfn} + endif + endif #------------------------------------------------------------ # Generate run script @@ -1144,7 +1189,9 @@ if (\${dobuild} == true) then endif endif if (\${dosubmit} == true) then - ./cice.submit | tee -a ../suite.jobs + set jobid = \`./cice.submit\` + echo "\$jobid" + echo "\$jobid \${ICE_TESTNAME} " >> ../suite.jobs else if (\${dorun} == true) then ./cice.test endif @@ -1159,6 +1206,10 @@ EOF echo "" endif +#set pd1 = `date -u "+%s%N"` +#@ pdd = ( $pd1 - $pd0 ) / 1000000 +#echo "tcxp case done $pdd" + # This is the foreach end for the testsuite end # This is the foreach end for the envnames @@ -1173,6 +1224,7 @@ if ( ${dosuite} == 1 ) then cat >> ${tsdir}/suite.submit << EOF0 set nonomatch && rm -f ciceexe.* && unset nonomatch +set nonomatch && rm -f ice_in_save* && unset nonomatch EOF0 @@ -1215,6 +1267,10 @@ endif #--------------------------------------------- +#set pd1 = `date -u "+%s%N"` +#@ pdd = ( $pd1 - $pd0 ) / 1000000 +#echo "tcxp done $pdd" + echo " " echo "${0} done" echo " " diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index 0ecc2ee5a..4b295b54d 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/cicecore/cicedynB/analysis/ice_history.F90 @@ -1624,6 +1624,8 @@ subroutine init_hist (dt) write(nu_diag,*) 'The following variables will be ', & 'written to the history tape: ' write(nu_diag,101) 'description','units','variable','frequency','x' + if (num_avail_hist_fields_tot == 0) & + write(nu_diag,*) '*** WARNING: NO HISTORY FIELDS WILL BE WRITTEN ***' do n=1,num_avail_hist_fields_tot if (avail_hist_fields(n)%vhistfreq_n /= 0) & write(nu_diag,100) avail_hist_fields(n)%vdesc, & @@ -1888,6 +1890,7 @@ subroutine accum_hist (dt) !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, & !$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt, & !$OMP worka,workb,worka3,Tinz4d,Sinz4d,Tsnz4d) + do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -1895,6 +1898,8 @@ subroutine accum_hist (dt) jlo = this_block%jlo jhi = this_block%jhi + if (allocated(a2D)) then + workb(:,:) = aice_init(:,:,iblk) ! if (f_example(1:1) /= 'x') & @@ -2879,6 +2884,10 @@ subroutine accum_hist (dt) call accum_hist_field(n_siforceintstry, iblk, worka(:,:), a2D) endif + endif ! if (allocated(a2D)) + + if (allocated(a3Dc)) then + ! 3D category fields if (f_aicen (1:1) /= 'x') & call accum_hist_field(n_aicen-n2D, iblk, ncat_hist, & @@ -2959,6 +2968,10 @@ subroutine accum_hist (dt) call accum_hist_field(n_siitdsnthick-n2D, iblk, ncat_hist, worka3(:,:,:), a3Dc) endif + endif ! if (allocated(a3Dc)) + + if (allocated(a4Di)) then + ! example for 3D field (x,y,z) ! if (f_field3dz (1:1) /= 'x') & ! call accum_hist_field(n_field3dz-n3Dccum, iblk, nzilyr, & @@ -2996,6 +3009,10 @@ subroutine accum_hist (dt) Sinz4d(:,:,1:nzilyr,1:ncat_hist), a4Di) endif + endif ! if (allocated(a3Dc)) + + if (allocated(a4Ds)) then + if (f_Tsnz (1:1) /= 'x') then Tsnz4d(:,:,:,:) = c0 do n = 1, ncat_hist @@ -3012,25 +3029,30 @@ subroutine accum_hist (dt) Tsnz4d(:,:,1:nzslyr,1:ncat_hist), a4Ds) endif - ! Calculate aggregate surface melt flux by summing category values - if (f_fmeltt_ai(1:1) /= 'x') then - do ns = 1, nstreams - if (n_fmeltt_ai(ns) /= 0) then - worka(:,:) = c0 - do j = jlo, jhi - do i = ilo, ihi - if (tmask(i,j,iblk)) then - do n=1,ncat_hist - worka(i,j) = worka(i,j) + a3Dc(i,j,n,n_fmelttn_ai(ns)-n2D,iblk) - enddo ! n - endif ! tmask - enddo ! i - enddo ! j - a2D(:,:,n_fmeltt_ai(ns),iblk) = worka(:,:) - endif - enddo - endif + endif ! if (allocated(a4Ds)) + + if (allocated(a3Dc) .and. allocated(a2D)) then + + ! Calculate aggregate surface melt flux by summing category values + if (f_fmeltt_ai(1:1) /= 'x') then + do ns = 1, nstreams + if (n_fmeltt_ai(ns) /= 0) then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + do n=1,ncat_hist + worka(i,j) = worka(i,j) + a3Dc(i,j,n,n_fmelttn_ai(ns)-n2D,iblk) + enddo ! n + endif ! tmask + enddo ! i + enddo ! j + a2D(:,:,n_fmeltt_ai(ns),iblk) = worka(:,:) + endif + enddo + endif + endif !--------------------------------------------------------------- ! accumulate other history output !--------------------------------------------------------------- diff --git a/cicecore/cicedynB/analysis/ice_history_bgc.F90 b/cicecore/cicedynB/analysis/ice_history_bgc.F90 index 67b23904e..fdb8c4393 100644 --- a/cicecore/cicedynB/analysis/ice_history_bgc.F90 +++ b/cicecore/cicedynB/analysis/ice_history_bgc.F90 @@ -778,9 +778,9 @@ subroutine init_hist_bgc_2D ! 2D variables - if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then + if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then - do ns = 1, nstreams + do ns = 1, nstreams if (histfreq(ns) /= 'x') then if (f_iso(1:1) /= 'x') then @@ -1782,9 +1782,9 @@ subroutine init_hist_bgc_2D ns, f_hbri) endif ! histfreq(ns) /= 'x' - enddo ! nstreams + enddo ! nstreams - endif ! tr_aero, etc + endif ! tr_aero, etc end subroutine init_hist_bgc_2D @@ -1841,7 +1841,7 @@ subroutine init_hist_bgc_3Db if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (z_tracers .or. solve_zsal) then + if (z_tracers .or. solve_zsal) then do ns = 1, nstreams if (histfreq(ns) /= 'x') then @@ -1880,7 +1880,7 @@ subroutine init_hist_bgc_3Db endif ! histfreq(ns) /= 'x' enddo ! ns - endif ! z_tracers or solve_zsal + endif ! z_tracers or solve_zsal end subroutine init_hist_bgc_3Db @@ -2017,10 +2017,10 @@ subroutine accum_hist_bgc (iblk) ! increment field !--------------------------------------------------------------- - ! 2d bgc fields - if (allocated(a2D)) then + ! 2d bgc fields + if (allocated(a2D)) then - if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then + if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then ! zsalinity if (f_fzsal (1:1) /= 'x') & @@ -2082,7 +2082,7 @@ subroutine accum_hist_bgc (iblk) enddo endif - if (skl_bgc) then + if (skl_bgc) then ! skeletal layer bgc @@ -2159,7 +2159,7 @@ subroutine accum_hist_bgc (iblk) call accum_hist_field(n_bgc_DMS, iblk, & sk_l*trcr(:,:,nt_bgc_DMS, iblk), a2D) - endif !skl_bgc + endif !skl_bgc ! skeletal layer and vertical bgc @@ -2354,7 +2354,7 @@ subroutine accum_hist_bgc (iblk) ! vertical biogeochemistry - if (z_tracers) then + if (z_tracers) then if (f_fzaero(1:1)/= 'x') then do n=1,n_zaero @@ -2634,30 +2634,30 @@ subroutine accum_hist_bgc (iblk) call accum_hist_field(n_PONfrac, iblk, & trcr(:,:,nt_zbgc_frac - 1 + nlt_bgc_PON, iblk), a2D) - endif ! z_tracers + endif ! z_tracers ! brine if (f_hbri (1:1) /= 'x') & call accum_hist_field(n_hbri, iblk, & hbri(:,:,iblk), a2D) - endif ! 2d bgc tracers, tr_aero, tr_brine, solve_zsal, skl_bgc - endif ! allocated(a2D) + endif ! 2d bgc tracers, tr_aero, tr_brine, solve_zsal, skl_bgc + endif ! allocated(a2D) ! 3D category fields - if (allocated(a3Dc)) then - if (tr_brine) then + if (allocated(a3Dc)) then + if (tr_brine) then ! 3Dc bgc category fields if (f_fbri (1:1) /= 'x') & call accum_hist_field(n_fbri-n2D, iblk, ncat_hist, & trcrn(:,:,nt_fbri,1:ncat_hist,iblk), a3Dc) - endif - endif ! allocated(a3Dc) + endif + endif ! allocated(a3Dc) - if (allocated(a3Db)) then - if (z_tracers .or. solve_zsal) then + if (allocated(a3Db)) then + if (z_tracers .or. solve_zsal) then ! 3Db category fields if (f_bTin (1:1) /= 'x') then @@ -2763,11 +2763,11 @@ subroutine accum_hist_bgc (iblk) workz(:,:,1:nzblyr), a3Db) endif - endif ! 3Db fields - endif ! allocated(a3Db) + endif ! 3Db fields + endif ! allocated(a3Db) - if (allocated(a3Da)) then - if (z_tracers) then + if (allocated(a3Da)) then + if (z_tracers) then ! 3Da category fields if (f_zaero (1:1) /= 'x') then @@ -3223,11 +3223,11 @@ subroutine init_hist_bgc_3Da if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - ! snow+bio grid - - if (z_tracers) then + ! snow+bio grid + + if (z_tracers) then - do ns = 1, nstreams + do ns = 1, nstreams if (histfreq(ns) /= 'x') then !---------------------------------------------------------------------------- diff --git a/cicecore/cicedynB/analysis/ice_history_drag.F90 b/cicecore/cicedynB/analysis/ice_history_drag.F90 index d2694fc9a..31a92158b 100644 --- a/cicecore/cicedynB/analysis/ice_history_drag.F90 +++ b/cicecore/cicedynB/analysis/ice_history_drag.F90 @@ -263,6 +263,8 @@ subroutine accum_hist_drag (iblk) ! 2D fields + if (allocated(a2D)) then + if (f_Cdn_atm (1:1) /= 'x') & call accum_hist_field(n_Cdn_atm, iblk, Cdn_atm(:,:,iblk), a2D) if (f_Cdn_ocn (1:1) /= 'x') & @@ -294,6 +296,7 @@ subroutine accum_hist_drag (iblk) iblk, Cdn_ocn_skin(:,:,iblk), a2D) end if + endif ! if(allocated(a2D)) endif ! formdrag end subroutine accum_hist_drag diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 9c52bb888..83374d4dd 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -326,11 +326,18 @@ subroutine eap (dt) do j = 1, ny_block do i = 1, nx_block if (icetmask(i,j,iblk)==0) then + if (tmask(i,j,iblk)) then ! structure tensor - a11_1(i,j,iblk) = p5 - a11_2(i,j,iblk) = p5 - a11_3(i,j,iblk) = p5 - a11_4(i,j,iblk) = p5 + a11_1(i,j,iblk) = p5 + a11_2(i,j,iblk) = p5 + a11_3(i,j,iblk) = p5 + a11_4(i,j,iblk) = p5 + else + a11_1(i,j,iblk) = c0 + a11_2(i,j,iblk) = c0 + a11_3(i,j,iblk) = c0 + a11_4(i,j,iblk) = c0 + endif a12_1(i,j,iblk) = c0 a12_2(i,j,iblk) = c0 a12_3(i,j,iblk) = c0 @@ -1932,58 +1939,38 @@ subroutine stepa (nx_block, ny_block, & j = indxtj(ij) ! ne - call calc_ffrac(1, stressp_1(i,j), stressm_1(i,j), & - stress12_1(i,j), & - a11_1(i,j), & - mresult11) - - call calc_ffrac(2, stressp_1(i,j), stressm_1(i,j), & - stress12_1(i,j), & - a12_1(i,j), & - mresult12) + call calc_ffrac(stressp_1(i,j), stressm_1(i,j), & + stress12_1(i,j), & + a11_1(i,j), a12_1(i,j), & + mresult11, mresult12) a11_1(i,j) = (a11_1(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit a12_1(i,j) = (a12_1(i,j)*dtei - mresult12) * dteikth ! implicit ! nw - call calc_ffrac(1, stressp_2(i,j), stressm_2(i,j), & - stress12_2(i,j), & - a11_2(i,j), & - mresult11) - - call calc_ffrac(2, stressp_2(i,j), stressm_2(i,j), & - stress12_2(i,j), & - a12_2(i,j), & - mresult12) + call calc_ffrac(stressp_2(i,j), stressm_2(i,j), & + stress12_2(i,j), & + a11_2(i,j), a12_2(i,j), & + mresult11, mresult12) a11_2(i,j) = (a11_2(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit a12_2(i,j) = (a12_2(i,j)*dtei - mresult12) * dteikth ! implicit ! sw - call calc_ffrac(1, stressp_3(i,j), stressm_3(i,j), & - stress12_3(i,j), & - a11_3(i,j), & - mresult11) - - call calc_ffrac(2, stressp_3(i,j), stressm_3(i,j), & - stress12_3(i,j), & - a12_3(i,j), & - mresult12) + call calc_ffrac(stressp_3(i,j), stressm_3(i,j), & + stress12_3(i,j), & + a11_3(i,j), a12_3(i,j), & + mresult11, mresult12) a11_3(i,j) = (a11_3(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit a12_3(i,j) = (a12_3(i,j)*dtei - mresult12) * dteikth ! implicit ! se - call calc_ffrac(1, stressp_4(i,j), stressm_4(i,j), & - stress12_4(i,j), & - a11_4(i,j), & - mresult11) - - call calc_ffrac(2, stressp_4(i,j), stressm_4(i,j), & - stress12_4(i,j), & - a12_4(i,j), & - mresult12) + call calc_ffrac(stressp_4(i,j), stressm_4(i,j), & + stress12_4(i,j), & + a11_4(i,j), a12_4(i,j), & + mresult11, mresult12) a11_4(i,j) = (a11_4(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit a12_4(i,j) = (a12_4(i,j)*dtei - mresult12) * dteikth ! implicit @@ -2003,19 +1990,17 @@ end subroutine stepa ! the ice floe re-orientation due to fracture ! Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2 - subroutine calc_ffrac (blockno, stressp, stressm, & - stress12, & - a1x, & - mresult) - integer(kind=int_kind), intent(in) :: & - blockno + subroutine calc_ffrac (stressp, stressm, & + stress12, & + a1x, a2x, & + mresult1, mresult2) real (kind=dbl_kind), intent(in) :: & - stressp, stressm, stress12, a1x + stressp, stressm, stress12, a1x, a2x real (kind=dbl_kind), intent(out) :: & - mresult + mresult1, mresult2 ! local variables @@ -2035,11 +2020,12 @@ subroutine calc_ffrac (blockno, stressp, stressm, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - sigma11 = p5*(stressp+stressm) - sigma12 = stress12 - sigma22 = p5*(stressp-stressm) + sigma11 = p5*(stressp+stressm) + sigma12 = stress12 + sigma22 = p5*(stressp-stressm) - if ((sigma11-sigma22) == c0) then +! if ((sigma11-sigma22) == c0) then sigma11-sigma22 == 0 => stressn ==0 + if (stressm == c0) then gamma = p5*(pih) else gamma = p5*atan2((c2*sigma12),(sigma11-sigma22)) @@ -2061,24 +2047,27 @@ subroutine calc_ffrac (blockno, stressp, stressm, & ! Pure divergence if ((sigma_1 >= c0).and.(sigma_2 >= c0)) then - mresult = c0 + mresult1 = c0 + mresult2 = c0 ! Unconfined compression: cracking of blocks not along the axial splitting direction ! which leads to the loss of their shape, so we again model it through diffusion elseif ((sigma_1 >= c0).and.(sigma_2 < c0)) then - if (blockno == 1) mresult = kfrac * (a1x - Q12Q12) - if (blockno == 2) mresult = kfrac * (a1x + Q11Q12) + mresult1 = kfrac * (a1x - Q12Q12) + mresult2 = kfrac * (a2x + Q11Q12) ! Shear faulting elseif (sigma_2 == c0) then - mresult = c0 + mresult1 = c0 + mresult2 = c0 elseif ((sigma_1 <= c0).and.(sigma_1/sigma_2 <= threshold)) then - if (blockno == 1) mresult = kfrac * (a1x - Q12Q12) - if (blockno == 2) mresult = kfrac * (a1x + Q11Q12) + mresult1 = kfrac * (a1x - Q12Q12) + mresult2 = kfrac * (a2x + Q11Q12) ! Horizontal spalling - else - mresult = c0 + else + mresult1 = c0 + mresult2 = c0 endif end subroutine calc_ffrac diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 276c8bb79..8f3fc4910 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -6,23 +6,25 @@ ! See: ! ! Hunke, E. C., and J. K. Dukowicz (1997). An elastic-viscous-plastic model -! for sea ice dynamics. {\em J. Phys. Oceanogr.}, {\bf 27}, 1849--1867. +! for sea ice dynamics. J. Phys. Oceanogr., 27, 1849-1867. ! ! Hunke, E. C. (2001). Viscous-Plastic Sea Ice Dynamics with the EVP Model: -! Linearization Issues. {\em Journal of Computational Physics}, {\bf 170}, -! 18--38. +! Linearization Issues. J. Comput. Phys., 170, 18-38. ! ! Hunke, E. C., and J. K. Dukowicz (2002). The Elastic-Viscous-Plastic ! Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates -! on a Sphere---Incorporation of Metric Terms. {\em Monthly Weather Review}, -! {\bf 130}, 1848--1865. +! on a Sphere - Incorporation of Metric Terms. Mon. Weather Rev., +! 130, 1848-1865. ! ! Hunke, E. C., and J. K. Dukowicz (2003). The sea ice momentum ! equation in the free drift regime. Los Alamos Tech. Rep. LA-UR-03-2219. ! -! Bouillon, S., T. Fichefet, V. Legat and G. Madec (submitted 2013). The -! revised elastic-viscous-plastic method. Ocean Modelling. +! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys. +! Oceanogr., 9, 817-846. ! +! Bouillon, S., T. Fichefet, V. Legat and G. Madec (2013). The +! elastic-viscous-plastic method revisited. Ocean Model., 71, 2-12. +! ! author: Elizabeth C. Hunke, LANL ! ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb (LANL) @@ -590,8 +592,8 @@ subroutine stress (nx_block, ny_block, & rdg_conv, rdg_shear, & str ) - use ice_dyn_shared, only: strain_rates, deformations - + use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure + integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ksub , & ! subcycling step @@ -640,9 +642,10 @@ subroutine stress (nx_block, ny_block, & tensionne, tensionnw, tensionse, tensionsw, & ! tension shearne, shearnw, shearse, shearsw , & ! shearing Deltane, Deltanw, Deltase, Deltasw , & ! Delt + zetax2ne, zetax2nw, zetax2se, zetax2sw , & ! 2 x zeta (visc coeff) + etax2ne, etax2nw, etax2se, etax2sw , & ! 2 x eta (visc coeff) + rep_prsne, rep_prsnw, rep_prsse, rep_prssw, & ! replacement pressure ! puny , & ! puny - c0ne, c0nw, c0se, c0sw , & ! useful combinations - c1ne, c1nw, c1se, c1sw , & ssigpn, ssigps, ssigpe, ssigpw , & ssigmn, ssigms, ssigme, ssigmw , & ssig12n, ssig12s, ssig12e, ssig12w , & @@ -653,6 +656,8 @@ subroutine stress (nx_block, ny_block, & str12ew, str12we, str12ns, str12sn , & strp_tmp, strm_tmp, tmp + logical :: capping ! of the viscous coef + character(len=*), parameter :: subname = '(stress)' !----------------------------------------------------------------- @@ -660,6 +665,7 @@ subroutine stress (nx_block, ny_block, & !----------------------------------------------------------------- str(:,:,:) = c0 + capping = .true. ! could be later included in ice_in do ij = 1, icellt i = indxti(ij) @@ -669,6 +675,7 @@ subroutine stress (nx_block, ny_block, & ! strain rates ! NOTE these are actually strain rates * area (m^2/s) !----------------------------------------------------------------- + call strain_rates (nx_block, ny_block, & i, j, & uvel, vvel, & @@ -685,46 +692,53 @@ subroutine stress (nx_block, ny_block, & Deltase, Deltasw ) !----------------------------------------------------------------- - ! strength/Delta ! kg/s + ! viscous coefficients and replacement pressure !----------------------------------------------------------------- - c0ne = strength(i,j)/max(Deltane,tinyarea(i,j)) - c0nw = strength(i,j)/max(Deltanw,tinyarea(i,j)) - c0sw = strength(i,j)/max(Deltasw,tinyarea(i,j)) - c0se = strength(i,j)/max(Deltase,tinyarea(i,j)) - - c1ne = c0ne*arlx1i - c1nw = c0nw*arlx1i - c1sw = c0sw*arlx1i - c1se = c0se*arlx1i - - c0ne = c1ne*ecci - c0nw = c1nw*ecci - c0sw = c1sw*ecci - c0se = c1se*ecci - + + call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),& + Deltane, Deltanw, & + Deltasw, Deltase, & + zetax2ne, zetax2nw, & + zetax2sw, zetax2se, & + etax2ne, etax2nw, & + etax2sw, etax2se, & + rep_prsne, rep_prsnw, & + rep_prssw, rep_prsse, & + capping) + !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + c1ne*(divune*(c1+Ktens) - Deltane*(c1-Ktens))) & - * denom1 - stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + c1nw*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens))) & - * denom1 - stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + c1sw*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens))) & - * denom1 - stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + c1se*(divuse*(c1+Ktens) - Deltase*(c1-Ktens))) & - * denom1 - - stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + c0ne*tensionne*(c1+Ktens)) * denom1 - stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + c0nw*tensionnw*(c1+Ktens)) * denom1 - stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + c0sw*tensionsw*(c1+Ktens)) * denom1 - stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + c0se*tensionse*(c1+Ktens)) * denom1 - - stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + c0ne*shearne*p5*(c1+Ktens)) * denom1 - stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + c0nw*shearnw*p5*(c1+Ktens)) * denom1 - stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + c0sw*shearsw*p5*(c1+Ktens)) * denom1 - stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + c0se*shearse*p5*(c1+Ktens)) * denom1 + ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code + + stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + & + arlx1i*(zetax2ne*divune - rep_prsne)) * denom1 + stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + & + arlx1i*(zetax2nw*divunw - rep_prsnw)) * denom1 + stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + & + arlx1i*(zetax2sw*divusw - rep_prssw)) * denom1 + stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + & + arlx1i*(zetax2se*divuse - rep_prsse)) * denom1 + + stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + & + arlx1i*etax2ne*tensionne) * denom1 + stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + & + arlx1i*etax2nw*tensionnw) * denom1 + stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + & + arlx1i*etax2sw*tensionsw) * denom1 + stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + & + arlx1i*etax2se*tensionse) * denom1 + + stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + & + arlx1i*p5*etax2ne*shearne) * denom1 + stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + & + arlx1i*p5*etax2nw*shearnw) * denom1 + stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + & + arlx1i*p5*etax2sw*shearsw) * denom1 + stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + & + arlx1i*p5*etax2se*shearse) * denom1 !----------------------------------------------------------------- ! Eliminate underflows. diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index bb65f122c..23251b2d1 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -26,6 +26,7 @@ module ice_dyn_shared dyn_prep1, dyn_prep2, dyn_finish, & seabed_stress_factor_LKD, seabed_stress_factor_prob, & alloc_dyn_shared, deformations, strain_rates, & + viscous_coeffs_and_rep_pressure, & stack_velocity_field, unstack_velocity_field ! namelist parameters @@ -62,8 +63,11 @@ 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 + e_yieldcurve, & ! VP aspect ratio of elliptical yield curve + e_plasticpot, & ! VP aspect ratio of elliptical plastic potential + epp2i , & ! 1/(e_plasticpot)^2 + e_factor , & ! (e_yieldcurve)^2/(e_plasticpot)^4 + ecci , & ! temporary for 1d evp dtei , & ! 1/dte, where dte is subcycling timestep (1/s) ! dte2T , & ! dte/2T denom1 ! constants for stress equation @@ -219,7 +223,6 @@ subroutine set_evp_parameters (dt) !real (kind=dbl_kind) :: & !dte , & ! subcycling timestep for EVP dynamics, s - !ecc , & ! (ratio of major to minor ellipse axes)^2 !tdamp2 ! 2*(wave damping time scale T) character(len=*), parameter :: subname = '(set_evp_parameters)' @@ -229,10 +232,10 @@ subroutine set_evp_parameters (dt) !dtei = c1/dte ! 1/s dtei = real(ndte,kind=dbl_kind)/dt - ! major/minor axis length ratio, squared - !ecc = e_ratio**2 - !ecci = c1/ecc ! 1/ecc - ecci = c1/e_ratio**2 ! 1/ecc + ! variables for elliptical yield curve and plastic potential + epp2i = c1/e_plasticpot**2 + e_factor = e_yieldcurve**2 / e_plasticpot**4 + ecci = c1/e_yieldcurve**2 ! temporary for 1d evp ! constants for stress equation !tdamp2 = c2*eyc*dt ! s @@ -864,7 +867,7 @@ end subroutine dyn_finish ! ! Lemieux, J. F., F. Dupont, P. Blain, F. Roy, G.C. Smith, G.M. Flato (2016). ! Improving the simulation of landfast ice by combining tensile strength and a -! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121. +! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121, 7354-7368. ! ! author: JF Lemieux, Philippe Blain (ECCC) ! @@ -1351,13 +1354,93 @@ subroutine strain_rates (nx_block, ny_block, & - cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j ) ! Delta (in the denominator of zeta, eta) - Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2)) - Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2)) - Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2)) - Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2)) + Deltane = sqrt(divune**2 + e_factor*(tensionne**2 + shearne**2)) + Deltanw = sqrt(divunw**2 + e_factor*(tensionnw**2 + shearnw**2)) + Deltasw = sqrt(divusw**2 + e_factor*(tensionsw**2 + shearsw**2)) + Deltase = sqrt(divuse**2 + e_factor*(tensionse**2 + shearse**2)) end subroutine strain_rates + !======================================================================= + ! Computes viscous coefficients and replacement pressure for stress + ! calculations. Note that tensile strength is included here. + ! + ! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys. + ! Oceanogr., 9, 817-846. + ! + ! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by + ! adding tensile strength. J. Phys. Oceanogr. 40, 185-198. + ! + ! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice + ! by combining tensile strength and a parameterization for grounded ridges. + ! J. Geophys. Res. Oceans, 121, 7354-7368. + + subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & + Deltane, Deltanw, & + Deltasw, Deltase, & + zetax2ne, zetax2nw, & + zetax2sw, zetax2se, & + etax2ne, etax2nw, & + etax2sw, etax2se, & + rep_prsne, rep_prsnw,& + rep_prssw, rep_prsse,& + capping) + + real (kind=dbl_kind), intent(in):: & + strength, tinyarea ! at the t-point + + real (kind=dbl_kind), intent(in):: & + Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner + + logical , intent(in):: capping + + real (kind=dbl_kind), intent(out):: & + zetax2ne, zetax2nw, zetax2sw, zetax2se, & ! zetax2 at each corner + etax2ne, etax2nw, etax2sw, etax2se, & ! etax2 at each corner + rep_prsne, rep_prsnw, rep_prssw, rep_prsse ! replacement pressure + + ! local variables + real (kind=dbl_kind) :: & + tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse + + ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code + +! if (trim(yield_curve) == 'ellipse') then + + if (capping) then + tmpcalcne = strength/max(Deltane,tinyarea) + tmpcalcnw = strength/max(Deltanw,tinyarea) + tmpcalcsw = strength/max(Deltasw,tinyarea) + tmpcalcse = strength/max(Deltase,tinyarea) + else + tmpcalcne = strength/(Deltane + tinyarea) + tmpcalcnw = strength/(Deltanw + tinyarea) + tmpcalcsw = strength/(Deltasw + tinyarea) + tmpcalcse = strength/(Deltase + tinyarea) + endif + + zetax2ne = (c1+Ktens)*tmpcalcne ! northeast + rep_prsne = (c1-Ktens)*tmpcalcne*Deltane + etax2ne = epp2i*zetax2ne + + zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest + rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw + etax2nw = epp2i*zetax2nw + + zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest + rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw + etax2sw = epp2i*zetax2sw + + zetax2se = (c1+Ktens)*tmpcalcse ! southeast + rep_prsse = (c1-Ktens)*tmpcalcse*Deltase + etax2se = epp2i*zetax2se + +! else + +! endif + + end subroutine viscous_coeffs_and_rep_pressure + !======================================================================= ! Load velocity components into array for boundary updates diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 860865dba..2f1285084 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -46,7 +46,7 @@ module ice_dyn_vp use ice_domain, only: nblocks, distrb_info use ice_domain_size, only: max_blocks use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, & - ecci, cosw, sinw, fcor_blk, uvel_init, vvel_init, & + cosw, sinw, fcor_blk, uvel_init, vvel_init, & seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & seabed_stress, Ktens, stack_velocity_field, unstack_velocity_field use ice_fileunits, only: nu_diag @@ -234,8 +234,10 @@ subroutine implicit_solver (dt) umassdti ! mass of U-cell/dte (kg/m^2 s) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4):: & - zetaD ! zetaD = 2zeta (viscous coeff) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + logical (kind=log_kind) :: calc_strair integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: & @@ -488,8 +490,9 @@ subroutine implicit_solver (dt) bxfix , byfix , & umassdti, sol , & fpresx , fpresy, & - zetaD , Cb , & - halo_info_mask) + zetax2 , etax2 , & + rep_prs , & + Cb, halo_info_mask) !----------------------------------------------------------------- ! End of nonlinear iteration !----------------------------------------------------------------- @@ -510,7 +513,8 @@ subroutine implicit_solver (dt) dxt (:,:,iblk), dyt (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & cxm (:,:,iblk), cym (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & + rep_prs (:,:,iblk,:), & stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & @@ -671,8 +675,9 @@ subroutine anderson_solver (icellt , icellu, & bxfix , byfix , & umassdti, sol , & fpresx , fpresy, & - zetaD , Cb , & - halo_info_mask) + zetax2 , etax2 , & + rep_prs , & + Cb, halo_info_mask) use ice_arrays_column, only: Cdn_ocn use ice_blocks, only: nx_block, ny_block @@ -708,8 +713,10 @@ subroutine anderson_solver (icellt , icellu, & umassdti ! mass of U-cell/dte (kg/m^2 s) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(out) :: & - zetaD ! zetaD = 2zeta (viscous coeff) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + type (ice_halo), intent(in) :: & halo_info_mask ! ghost cell update info for masked halo @@ -805,7 +812,7 @@ subroutine anderson_solver (icellt , icellu, & do it_nl = 0, maxits_nonlin ! nonlinear iteration loop ! Compute quantities needed for computing PDE residual and g(x) (fixed point map) !----------------------------------------------------------------- - ! Calc zetaD, dPr/dx, dPr/dy, Cb and vrel = f(uprev_k, vprev_k) + ! Calc zetax2, etax2, dPr/dx, dPr/dy, Cb and vrel = f(uprev_k, vprev_k) !----------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks @@ -828,9 +835,9 @@ subroutine anderson_solver (icellt , icellu, & dxhy (:,:,iblk), dyhx (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & cxm (:,:,iblk), cym (:,:,iblk), & - tinyarea (:,:,iblk), & - strength (:,:,iblk), zetaD (:,:,iblk,:), & - stress_Pr (:,:,:)) + tinyarea (:,:,iblk), strength (:,:,iblk),& + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:),& + rep_prs(:,:,iblk,:), stress_Pr (:,:,:)) call calc_vrel_Cb (nx_block , ny_block , & icellu (iblk), Cdn_ocn (:,:,iblk), & @@ -861,7 +868,7 @@ subroutine anderson_solver (icellt , icellu, & cxm (:,:,iblk) , cym (:,:,iblk), & uprev_k (:,:,iblk) , vprev_k (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & Au (:,:,iblk) , Av (:,:,iblk)) @@ -908,14 +915,15 @@ subroutine anderson_solver (icellt , icellu, & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks ! first compute diagonal contributions due to rheology term - call formDiag_step1 (nx_block , ny_block , & - icellu (iblk) , & - indxui (:,iblk) , indxuj(:,iblk), & - dxt (:,:,iblk) , dyt (:,:,iblk), & - dxhy (:,:,iblk) , dyhx(:,:,iblk), & - cxp (:,:,iblk) , cyp (:,:,iblk), & - cxm (:,:,iblk) , cym (:,:,iblk), & - zetaD (:,:,iblk,:), diag_rheo(:,:,:)) + call formDiag_step1 (nx_block , ny_block , & + icellu (iblk) , & + indxui (:,iblk) , indxuj(:,iblk), & + dxt (:,:,iblk) , dyt (:,:,iblk), & + dxhy (:,:,iblk) , dyhx(:,:,iblk), & + cxp (:,:,iblk) , cyp (:,:,iblk), & + cxm (:,:,iblk) , cym (:,:,iblk), & + zetax2 (:,:,iblk,:), etax2(:,:,iblk,:), & + diag_rheo(:,:,:)) ! second compute the full diagonal call formDiag_step2 (nx_block , ny_block , & icellu (iblk), & @@ -929,7 +937,7 @@ subroutine anderson_solver (icellt , icellu, & endif ! FGMRES linear solver - call fgmres (zetaD , & + call fgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & halo_info_mask, & @@ -1119,7 +1127,7 @@ end subroutine anderson_solver !======================================================================= -! Computes the viscous coefficients (in fact zetaD=2*zeta) and dPr/dx, dPr/dy +! Computes the viscous coefficients and dPr/dx, dPr/dy subroutine calc_zeta_dPr (nx_block, ny_block, & icellt , & @@ -1129,11 +1137,11 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & dxhy , dyhx , & cxp , cyp , & cxm , cym , & - tinyarea, & - strength, zetaD , & - stPr) + tinyarea, strength, & + zetax2 , etax2 , & + rep_prs , stPr) - use ice_dyn_shared, only: strain_rates + use ice_dyn_shared, only: strain_rates, viscous_coeffs_and_rep_pressure integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1158,8 +1166,10 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & tinyarea ! min_strain_rate*tarea real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(out) :: & - zetaD ! 2*zeta - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: & stPr ! stress combinations from replacement pressure @@ -1186,9 +1196,11 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & capping = .false. - ! Initialize stPr and zetaD to zero (for cells where icetmask is false) - stPr = c0 - zetaD = c0 + ! Initialize stPr, zetax2 and etax2 to zero + ! (for cells where icetmask is false) + stPr = c0 + zetax2 = c0 + etax2 = c0 do ij = 1, icellt i = indxti(ij) @@ -1213,27 +1225,30 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & Deltane , Deltanw , & Deltase , Deltasw) - if (capping) then - zetaD(i,j,1) = strength(i,j)/max(Deltane,tinyarea(i,j)) - zetaD(i,j,2) = strength(i,j)/max(Deltanw,tinyarea(i,j)) - zetaD(i,j,3) = strength(i,j)/max(Deltasw,tinyarea(i,j)) - zetaD(i,j,4) = strength(i,j)/max(Deltase,tinyarea(i,j)) - else - zetaD(i,j,1) = strength(i,j)/(Deltane + tinyarea(i,j)) - zetaD(i,j,2) = strength(i,j)/(Deltanw + tinyarea(i,j)) - zetaD(i,j,3) = strength(i,j)/(Deltasw + tinyarea(i,j)) - zetaD(i,j,4) = strength(i,j)/(Deltase + tinyarea(i,j)) - endif + !----------------------------------------------------------------- + ! viscous coefficients and replacement pressure + !----------------------------------------------------------------- + + call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j), & + Deltane, Deltanw, & + Deltasw, Deltase, & + zetax2(i,j,1), zetax2(i,j,2), & + zetax2(i,j,3), zetax2(i,j,4), & + etax2(i,j,1), etax2(i,j,2), & + etax2(i,j,3), etax2(i,j,4), & + rep_prs(i,j,1), rep_prs(i,j,2), & + rep_prs(i,j,3), rep_prs(i,j,4), & + capping) !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1 = -zetaD(i,j,1)*(Deltane*(c1-Ktens)) - stressp_2 = -zetaD(i,j,2)*(Deltanw*(c1-Ktens)) - stressp_3 = -zetaD(i,j,3)*(Deltasw*(c1-Ktens)) - stressp_4 = -zetaD(i,j,4)*(Deltase*(c1-Ktens)) + stressp_1 = -rep_prs(i,j,1) + stressp_2 = -rep_prs(i,j,2) + stressp_3 = -rep_prs(i,j,3) + stressp_4 = -rep_prs(i,j,4) !----------------------------------------------------------------- ! combinations of the Pr related stresses for the momentum equation ! kg/s^2 @@ -1305,6 +1320,9 @@ end subroutine calc_zeta_dPr ! Computes the VP stresses (as diagnostic) +! Lemieux, J.-F., and Dupont, F. (2020), On the calculation of normalized +! viscous-plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769, + subroutine stress_vp (nx_block , ny_block , & icellt , & indxti , indxtj , & @@ -1312,7 +1330,8 @@ subroutine stress_vp (nx_block , ny_block , & dxt , dyt , & cxp , cyp , & cxm , cym , & - zetaD , & + zetax2 , etax2 , & + rep_prs , & stressp_1 , stressp_2 , & stressp_3 , stressp_4 , & stressm_1 , stressm_2 , & @@ -1341,8 +1360,10 @@ subroutine stress_vp (nx_block , ny_block , & cxm ! 0.5*HTN - 1.5*HTS real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 @@ -1388,21 +1409,21 @@ subroutine stress_vp (nx_block , ny_block , & ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - - stressp_1(i,j) = zetaD(i,j,1)*(divune*(c1+Ktens) - Deltane*(c1-Ktens)) - stressp_2(i,j) = zetaD(i,j,2)*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens)) - stressp_3(i,j) = zetaD(i,j,3)*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens)) - stressp_4(i,j) = zetaD(i,j,4)*(divuse*(c1+Ktens) - Deltase*(c1-Ktens)) - stressm_1(i,j) = zetaD(i,j,1)*tensionne*(c1+Ktens)*ecci - stressm_2(i,j) = zetaD(i,j,2)*tensionnw*(c1+Ktens)*ecci - stressm_3(i,j) = zetaD(i,j,3)*tensionsw*(c1+Ktens)*ecci - stressm_4(i,j) = zetaD(i,j,4)*tensionse*(c1+Ktens)*ecci + stressp_1(i,j) = zetax2(i,j,1)*divune - rep_prs(i,j,1) + stressp_2(i,j) = zetax2(i,j,2)*divunw - rep_prs(i,j,2) + stressp_3(i,j) = zetax2(i,j,3)*divusw - rep_prs(i,j,3) + stressp_4(i,j) = zetax2(i,j,4)*divuse - rep_prs(i,j,4) - stress12_1(i,j) = zetaD(i,j,1)*shearne*p5*(c1+Ktens)*ecci - stress12_2(i,j) = zetaD(i,j,2)*shearnw*p5*(c1+Ktens)*ecci - stress12_3(i,j) = zetaD(i,j,3)*shearsw*p5*(c1+Ktens)*ecci - stress12_4(i,j) = zetaD(i,j,4)*shearse*p5*(c1+Ktens)*ecci + stressm_1(i,j) = etax2(i,j,1)*tensionne + stressm_2(i,j) = etax2(i,j,2)*tensionnw + stressm_3(i,j) = etax2(i,j,3)*tensionsw + stressm_4(i,j) = etax2(i,j,4)*tensionse + + stress12_1(i,j) = etax2(i,j,1)*shearne*p5 + stress12_2(i,j) = etax2(i,j,2)*shearnw*p5 + stress12_3(i,j) = etax2(i,j,3)*shearsw*p5 + stress12_4(i,j) = etax2(i,j,4)*shearse*p5 enddo ! ij @@ -1534,7 +1555,7 @@ subroutine matvec (nx_block, ny_block, & cxm , cym , & uvel , vvel , & vrel , Cb , & - zetaD , & + zetax2 , etax2 , & umassdti, fm , & uarear , & Au , Av) @@ -1572,7 +1593,8 @@ subroutine matvec (nx_block, ny_block, & uarear ! 1/uarea real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & Au , & ! matvec, Fx = bx - Au (N/m^2) @@ -1647,20 +1669,20 @@ subroutine matvec (nx_block, ny_block, & ! NOTE: commented part of stressp is from the replacement pressure Pr !----------------------------------------------------------------- - stressp_1 = zetaD(i,j,1)*(divune*(c1+Ktens))! - Deltane*(c1-Ktens)) - stressp_2 = zetaD(i,j,2)*(divunw*(c1+Ktens))! - Deltanw*(c1-Ktens)) - stressp_3 = zetaD(i,j,3)*(divusw*(c1+Ktens))! - Deltasw*(c1-Ktens)) - stressp_4 = zetaD(i,j,4)*(divuse*(c1+Ktens))! - Deltase*(c1-Ktens)) + stressp_1 = zetax2(i,j,1)*divune! - Deltane*(c1-Ktens)) + stressp_2 = zetax2(i,j,2)*divunw! - Deltanw*(c1-Ktens)) + stressp_3 = zetax2(i,j,3)*divusw! - Deltasw*(c1-Ktens)) + stressp_4 = zetax2(i,j,4)*divuse! - Deltase*(c1-Ktens)) - stressm_1 = zetaD(i,j,1)*tensionne*(c1+Ktens)*ecci - stressm_2 = zetaD(i,j,2)*tensionnw*(c1+Ktens)*ecci - stressm_3 = zetaD(i,j,3)*tensionsw*(c1+Ktens)*ecci - stressm_4 = zetaD(i,j,4)*tensionse*(c1+Ktens)*ecci + stressm_1 = etax2(i,j,1)*tensionne + stressm_2 = etax2(i,j,2)*tensionnw + stressm_3 = etax2(i,j,3)*tensionsw + stressm_4 = etax2(i,j,4)*tensionse - stress12_1 = zetaD(i,j,1)*shearne*p5*(c1+Ktens)*ecci - stress12_2 = zetaD(i,j,2)*shearnw*p5*(c1+Ktens)*ecci - stress12_3 = zetaD(i,j,3)*shearsw*p5*(c1+Ktens)*ecci - stress12_4 = zetaD(i,j,4)*shearse*p5*(c1+Ktens)*ecci + stress12_1 = etax2(i,j,1)*shearne*p5 + stress12_2 = etax2(i,j,2)*shearnw*p5 + stress12_3 = etax2(i,j,3)*shearsw*p5 + stress12_4 = etax2(i,j,4)*shearse*p5 !----------------------------------------------------------------- ! combinations of the stresses for the momentum equation ! kg/s^2 @@ -1991,7 +2013,8 @@ subroutine formDiag_step1 (nx_block, ny_block, & dxhy , dyhx , & cxp , cyp , & cxm , cym , & - zetaD , Drheo) + zetax2 , etax2 , & + Drheo) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -2012,7 +2035,8 @@ subroutine formDiag_step1 (nx_block, ny_block, & cxm ! 0.5*HTN - 1.5*HTS real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: & Drheo ! intermediate value for diagonal components of matrix A associated @@ -2200,20 +2224,20 @@ subroutine formDiag_step1 (nx_block, ny_block, & ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1 = zetaD(i,j,1)*divune*(c1+Ktens) - stressp_2 = zetaD(i,j,2)*divunw*(c1+Ktens) - stressp_3 = zetaD(i,j,3)*divusw*(c1+Ktens) - stressp_4 = zetaD(i,j,4)*divuse*(c1+Ktens) + stressp_1 = zetax2(i,j,1)*divune + stressp_2 = zetax2(i,j,2)*divunw + stressp_3 = zetax2(i,j,3)*divusw + stressp_4 = zetax2(i,j,4)*divuse - stressm_1 = zetaD(i,j,1)*tensionne*(c1+Ktens)*ecci - stressm_2 = zetaD(i,j,2)*tensionnw*(c1+Ktens)*ecci - stressm_3 = zetaD(i,j,3)*tensionsw*(c1+Ktens)*ecci - stressm_4 = zetaD(i,j,4)*tensionse*(c1+Ktens)*ecci + stressm_1 = etax2(i,j,1)*tensionne + stressm_2 = etax2(i,j,2)*tensionnw + stressm_3 = etax2(i,j,3)*tensionsw + stressm_4 = etax2(i,j,4)*tensionse - stress12_1 = zetaD(i,j,1)*shearne*p5*(c1+Ktens)*ecci - stress12_2 = zetaD(i,j,2)*shearnw*p5*(c1+Ktens)*ecci - stress12_3 = zetaD(i,j,3)*shearsw*p5*(c1+Ktens)*ecci - stress12_4 = zetaD(i,j,4)*shearse*p5*(c1+Ktens)*ecci + stress12_1 = etax2(i,j,1)*shearne*p5 + stress12_2 = etax2(i,j,2)*shearnw*p5 + stress12_3 = etax2(i,j,3)*shearsw*p5 + stress12_4 = etax2(i,j,4)*shearse*p5 !----------------------------------------------------------------- ! combinations of the stresses for the momentum equation ! kg/s^2 @@ -2657,7 +2681,7 @@ end subroutine qr_delete ! ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC - subroutine fgmres (zetaD , & + subroutine fgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & halo_info_mask , & @@ -2673,8 +2697,9 @@ subroutine fgmres (zetaD , & use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -2784,7 +2809,7 @@ subroutine fgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & solx (:,:,iblk) , soly (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) @@ -2855,7 +2880,7 @@ subroutine fgmres (zetaD , & initer = initer + 1 nextit = initer + 1 ! precondition the current Arnoldi vector - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & arnoldi_basis_x(:,:,:,initer), & @@ -2891,7 +2916,7 @@ subroutine fgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & arnoldi_basis_x(:,:,iblk,nextit), & @@ -3057,7 +3082,7 @@ end subroutine fgmres ! ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC - subroutine pgmres (zetaD , & + subroutine pgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & bx , by , & @@ -3068,8 +3093,9 @@ subroutine pgmres (zetaD , & nbiter) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -3176,7 +3202,7 @@ subroutine pgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & solx (:,:,iblk) , soly (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) @@ -3248,7 +3274,7 @@ subroutine pgmres (zetaD , & nextit = initer + 1 ! precondition the current Arnoldi vector - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & arnoldi_basis_x(:,:,:,initer), & @@ -3272,7 +3298,7 @@ subroutine pgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & arnoldi_basis_x(:,:,iblk,nextit), & @@ -3385,7 +3411,7 @@ subroutine pgmres (zetaD , & end do ! Call preconditioner - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & workspace_x , workspace_y, & @@ -3451,7 +3477,7 @@ end subroutine pgmres ! ! authors: Philippe Blain, ECCC - subroutine precondition(zetaD , & + subroutine precondition(zetax2 , etax2, & Cb , vrel , & umassdti , & vx , vy , & @@ -3460,8 +3486,9 @@ subroutine precondition(zetaD , & wx , wy) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -3525,7 +3552,7 @@ subroutine precondition(zetaD , & tolerance = reltol_pgmres maxinner = dim_pgmres maxouter = maxits_pgmres - call pgmres (zetaD, & + call pgmres (zetax2, etax2 , & Cb , vrel , & umassdti , & vx , vy , & diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index 84bf1d461..cede58950 100755 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -118,7 +118,7 @@ module ice_forcing atm_data_format, & ! 'bin'=binary or 'nc'=netcdf ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf atm_data_type, & ! 'default', 'monthly', 'ncar', - ! 'LYq' or 'hadgem' or 'oned' or + ! 'hadgem' or 'oned' or ! 'JRA55_gx1' or 'JRA55_gx3' or 'JRA55_tx1' bgc_data_type, & ! 'default', 'clim' ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned', @@ -290,8 +290,10 @@ subroutine init_forcing_atmo ! default forcing values from init_flux_atm if (trim(atm_data_type) == 'ncar') then call NCAR_files(fyear) +#ifdef UNDEPRECATE_LYq elseif (trim(atm_data_type) == 'LYq') then call LY_files(fyear) +#endif elseif (trim(atm_data_type) == 'JRA55_gx1') then call JRA55_gx1_files(fyear) elseif (trim(atm_data_type) == 'JRA55_gx3') then @@ -606,8 +608,10 @@ subroutine get_forcing_atmo if (trim(atm_data_type) == 'ncar') then call ncar_data +#ifdef UNDEPRECATE_LYq elseif (trim(atm_data_type) == 'LYq') then call LY_data +#endif elseif (trim(atm_data_type) == 'JRA55_gx1') then call JRA55_data elseif (trim(atm_data_type) == 'JRA55_gx3') then @@ -1621,6 +1625,7 @@ subroutine prepare_forcing (nx_block, ny_block, & enddo enddo +#ifdef UNDEPRECATE_LYq elseif (trim(atm_data_type) == 'LYq') then ! precip is in mm/s @@ -1636,7 +1641,7 @@ subroutine prepare_forcing (nx_block, ny_block, & hm(i,j), flw(i,j)) enddo enddo - +#endif elseif (trim(atm_data_type) == 'oned') then ! rectangular grid ! precip is in kg/m^2/s @@ -2089,6 +2094,7 @@ subroutine ncar_data end subroutine ncar_data +#ifdef UNDEPRECATE_LYq !======================================================================= ! Large and Yeager forcing (AOMIP style) !======================================================================= @@ -2145,7 +2151,7 @@ subroutine LY_files (yr) endif ! master_task end subroutine LY_files - +#endif !======================================================================= subroutine JRA55_gx1_files(yr) @@ -2209,6 +2215,7 @@ subroutine JRA55_gx3_files(yr) endif end subroutine JRA55_gx3_files +#ifdef UNDEPRECATE_LYq !======================================================================= ! ! read Large and Yeager atmospheric data @@ -2432,7 +2439,7 @@ subroutine LY_data endif ! debug_forcing end subroutine LY_data - +#endif !======================================================================= subroutine JRA55_data diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 3d102217a..7485cbe23 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -105,9 +105,10 @@ subroutine input_data use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, & evp_algorithm, & seabed_stress, seabed_stress_method, & - k1, k2, alphab, threshold_hw, & - Ktens, e_ratio, coriolis, ssh_stress, & - kridge, brlx, arlx + k1, k2, alphab, threshold_hw, Ktens, & + e_yieldcurve, e_plasticpot, coriolis, & + ssh_stress, kridge, brlx, arlx + use ice_dyn_vp, only: maxits_nonlin, precond, dim_fgmres, dim_pgmres, maxits_fgmres, & maxits_pgmres, monitor_nonlin, monitor_fgmres, & monitor_pgmres, reltol_nonlin, reltol_fgmres, reltol_pgmres, & @@ -208,20 +209,20 @@ subroutine input_data namelist /dynamics_nml/ & kdyn, ndte, revised_evp, yield_curve, & - evp_algorithm, & + evp_algorithm, & brlx, arlx, ssh_stress, & advection, coriolis, kridge, ktransport, & kstrength, krdg_partic, krdg_redist, mu_rdg, & - e_ratio, Ktens, Cf, seabed_stress, & - k1, maxits_nonlin, precond, dim_fgmres, & + e_yieldcurve, e_plasticpot, Ktens, & + maxits_nonlin, precond, dim_fgmres, & dim_pgmres, maxits_fgmres, maxits_pgmres, monitor_nonlin, & monitor_fgmres, monitor_pgmres, reltol_nonlin, reltol_fgmres, & reltol_pgmres, algo_nonlin, dim_andacc, reltol_andacc, & damping_andacc, start_andacc, fpfunc_andacc, use_mean_vrel, & - ortho_type, & - k2, alphab, threshold_hw, & - seabed_stress_method, Pstar, Cstar - + ortho_type, seabed_stress, seabed_stress_method, & + k1, k2, alphab, threshold_hw, & + Cf, Pstar, Cstar + namelist /shortwave_nml/ & shortwave, albedo_type, & albicev, albicei, albsnowv, albsnowi, & @@ -367,7 +368,8 @@ subroutine input_data alphab = 20.0_dbl_kind ! alphab=Cb factor in Lemieux et al 2015 threshold_hw = 30.0_dbl_kind ! max water depth for grounding Ktens = 0.0_dbl_kind ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) - e_ratio = 2.0_dbl_kind ! VP ellipse aspect ratio + e_yieldcurve = 2.0_dbl_kind ! VP aspect ratio of elliptical yield curve + e_plasticpot = 2.0_dbl_kind ! VP aspect ratio of elliptical plastic potential maxits_nonlin = 4 ! max nb of iteration for nonlinear solver precond = 'pgmres' ! preconditioner for fgmres: 'ident' (identity), 'diag' (diagonal), 'pgmres' (Jacobi-preconditioned GMRES) dim_fgmres = 50 ! size of fgmres Krylov subspace @@ -448,7 +450,7 @@ subroutine input_data albsnowv = 0.98_dbl_kind ! cold snow albedo, visible albsnowi = 0.70_dbl_kind ! cold snow albedo, near IR ahmax = 0.3_dbl_kind ! thickness above which ice albedo is constant (m) - atmbndy = 'default' ! or 'constant' + atmbndy = 'similarity' ! Atm boundary layer: 'similarity', 'constant' or 'mixed' default_season = 'winter' ! default forcing data, if data is not read in fyear_init = 1900 ! first year of forcing cycle ycycle = 1 ! number of years in forcing cycle @@ -729,7 +731,8 @@ subroutine input_data call broadcast_scalar(alphab, master_task) call broadcast_scalar(threshold_hw, master_task) call broadcast_scalar(Ktens, master_task) - call broadcast_scalar(e_ratio, master_task) + call broadcast_scalar(e_yieldcurve, master_task) + call broadcast_scalar(e_plasticpot, master_task) call broadcast_scalar(advection, master_task) call broadcast_scalar(conserv_check, master_task) call broadcast_scalar(shortwave, master_task) @@ -955,6 +958,10 @@ subroutine input_data abort_list = trim(abort_list)//":1" endif + if (ktransport <= 0) then + advection = 'none' + endif + if (ktransport > 0 .and. advection /= 'remap' .and. advection /= 'upwind') then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: invalid advection=',trim(advection) abort_list = trim(abort_list)//":3" @@ -1214,6 +1221,14 @@ subroutine input_data endif endif + if (trim(atmbndy) == 'default') then + if (my_task == master_task) then + write(nu_diag,*) subname//' WARNING: atmbndy = default is deprecated' + write(nu_diag,*) subname//' WARNING: setting atmbndy = similarity' + endif + atmbndy = 'similarity' + endif + if (formdrag) then if (trim(atmbndy) == 'constant') then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: formdrag=T and atmbndy=constant' @@ -1432,9 +1447,10 @@ subroutine input_data if (kdyn == 1 .or. kdyn == 3) then write(nu_diag,1030) ' yield_curve = ', trim(yield_curve) if (trim(yield_curve) == 'ellipse') & - write(nu_diag,1002) ' e_ratio = ', e_ratio, ' : aspect ratio of ellipse' + write(nu_diag,1002) ' e_yieldcurve = ', e_yieldcurve, ' : aspect ratio of yield curve' + write(nu_diag,1002) ' e_plasticpot = ', e_plasticpot, ' : aspect ratio of plastic potential' endif - + if (trim(coriolis) == 'latitude') then tmpstr2 = ' : latitude-dependent Coriolis parameter' elseif (trim(coriolis) == 'contant') then @@ -1455,9 +1471,6 @@ subroutine input_data endif write(nu_diag,1030) ' ssh_stress = ',trim(ssh_stress),trim(tmpstr2) - if (ktransport <= 0) then - advection = 'none' - endif if (trim(advection) == 'remap') then tmpstr2 = ' : linear remapping advection' elseif (trim(advection) == 'upwind') then @@ -1641,13 +1654,18 @@ subroutine input_data write(nu_diag,1010) ' rotate_wind = ', rotate_wind,' : rotate wind/stress to computational grid' write(nu_diag,1010) ' formdrag = ', formdrag,' : use form drag parameterization' write(nu_diag,1000) ' iceruf = ', iceruf, ' : ice surface roughness at atmosphere interface (m)' - if (trim(atmbndy) == 'default') then - tmpstr2 = ' : stability-based boundary layer' + if (trim(atmbndy) == 'constant') then + tmpstr2 = ' : constant-based boundary layer' + elseif (trim(atmbndy) == 'similarity' .or. & + trim(atmbndy) == 'mixed') then write(nu_diag,1010) ' highfreq = ', highfreq,' : high-frequency atmospheric coupling' write(nu_diag,1020) ' natmiter = ', natmiter,' : number of atmo boundary layer iterations' write(nu_diag,1002) ' atmiter_conv = ', atmiter_conv,' : convergence criterion for ustar' - elseif (trim(atmbndy) == 'constant') then - tmpstr2 = ' : boundary layer uses bulk transfer coefficients' + if (trim(atmbndy) == 'similarity') then + tmpstr2 = ' : stability-based boundary layer' + else + tmpstr2 = ' : stability-based boundary layer for wind stress, constant-based for sensible+latent heat fluxes' + endif else tmpstr2 = ' : unknown value' endif diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index cfca994c3..8b69730b8 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -8,6 +8,7 @@ module CICE_InitMod use icepack_intfc, only: icepack_aggregate use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_init_snow use icepack_intfc, only: icepack_configure use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags @@ -83,7 +84,7 @@ subroutine cice_init2() use ice_dyn_vp , only: init_vp use ice_flux , only: init_coupler_flux, init_history_therm use ice_flux , only: init_history_dyn, init_flux_atm, init_flux_ocn - use ice_forcing , only: init_forcing_ocn + use ice_forcing , only: init_forcing_ocn, init_snowtable use ice_forcing_bgc , only: get_forcing_bgc, get_atm_bgc use ice_forcing_bgc , only: faero_default, faero_optics, alloc_forcing_bgc, fiso_default use ice_history , only: init_hist, accum_hist @@ -95,7 +96,8 @@ subroutine cice_init2() use ice_transport_driver , only: init_transport logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers - logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec + logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow + character(len=char_len) :: snw_aging_table character(len=*), parameter :: subname = '(cice_init2)' !---------------------------------------------------- @@ -145,7 +147,7 @@ subroutine cice_init2() call ice_HaloRestore_init ! restored boundary conditions call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & - wave_spec_out=wave_spec) + wave_spec_out=wave_spec, snw_aging_table_out=snw_aging_table) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(trim(subname), & file=__FILE__,line= __LINE__) @@ -158,7 +160,7 @@ subroutine cice_init2() call init_history_dyn ! initialize dynamic history variables call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) - call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_query_tracer_flags(tr_iso_out=tr_iso, tr_snow_out=tr_snow) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(trim(subname), & file=__FILE__,line= __LINE__) @@ -167,6 +169,17 @@ subroutine cice_init2() call faero_optics !initialize aerosol optical property tables end if + ! snow aging lookup table initialization + if (tr_snow) then ! advanced snow physics + call icepack_init_snow() + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + if (snw_aging_table(1:4) /= 'test') then + call init_snowtable() + endif + endif + ! Initialize shortwave components using swdn from previous timestep ! if restarting. These components will be scaled to current forcing ! in prep_radiation. @@ -199,12 +212,12 @@ subroutine init_restart() use ice_calendar, only: calendar use ice_constants, only: c0 use ice_domain, only: nblocks - use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd, nslyr use ice_dyn_eap, only: read_restart_eap use ice_dyn_shared, only: kdyn use ice_grid, only: tmask use ice_init, only: ice_ic - use ice_init_column, only: init_age, init_FY, init_lvl, & + use ice_init_column, only: init_age, init_FY, init_lvl, init_snowtracers, & init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd use ice_restart_column, only: restart_age, read_restart_age, & @@ -212,6 +225,7 @@ subroutine init_restart() restart_pond_cesm, read_restart_pond_cesm, & restart_pond_lvl, read_restart_pond_lvl, & restart_pond_topo, read_restart_pond_topo, & + restart_snow, read_restart_snow, & restart_fsd, read_restart_fsd, & restart_iso, read_restart_iso, & restart_aero, read_restart_aero, & @@ -226,12 +240,13 @@ subroutine init_restart() iblk ! block index logical(kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & - tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, tr_snow, & skl_bgc, z_tracers, solve_zsal integer(kind=int_kind) :: & ntrcr integer(kind=int_kind) :: & nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_smice, nt_smliq, nt_rhos, nt_rsnw, & nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice character(len=*), parameter :: subname = '(init_restart)' @@ -247,10 +262,12 @@ subroutine init_restart() call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & - tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + tr_snow_out=tr_snow, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & + nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw, & nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & @@ -347,6 +364,21 @@ subroutine init_restart() enddo ! iblk endif ! .not. restart_pond endif + ! snow redistribution/metamorphism + if (tr_snow) then + if (trim(runtype) == 'continue') restart_snow = .true. + if (restart_snow) then + call read_restart_snow + else + do iblk = 1, nblocks + call init_snowtracers(trcrn(:,:,nt_smice:nt_smice+nslyr-1,:,iblk), & + trcrn(:,:,nt_smliq:nt_smliq+nslyr-1,:,iblk), & + trcrn(:,:,nt_rhos :nt_rhos +nslyr-1,:,iblk), & + trcrn(:,:,nt_rsnw :nt_rsnw +nslyr-1,:,iblk)) + enddo ! iblk + endif + endif + ! floe size distribution if (tr_fsd) then if (trim(runtype) == 'continue') restart_fsd = .true. @@ -441,7 +473,6 @@ subroutine init_restart() call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - end subroutine init_restart !======================================================================= diff --git a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 index 81fa367c1..219777f6f 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 @@ -110,7 +110,7 @@ subroutine ice_step use ice_boundary, only: ice_HaloUpdate use ice_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep use ice_calendar, only: idate, msec - use ice_diagnostics, only: init_mass_diags, runtime_diags + use ice_diagnostics, only: init_mass_diags, runtime_diags, debug_model, debug_ice use ice_diagnostics_bgc, only: hbrine_diags, zsal_diags, bgc_diags use ice_domain, only: halo_info, nblocks use ice_dyn_eap, only: write_restart_eap @@ -123,12 +123,13 @@ subroutine ice_step use ice_restart_column, only: write_restart_age, write_restart_FY, & write_restart_lvl, write_restart_pond_cesm, write_restart_pond_lvl, & write_restart_pond_topo, write_restart_aero, write_restart_fsd, & - write_restart_iso, write_restart_bgc, write_restart_hbrine + write_restart_iso, write_restart_bgc, write_restart_hbrine, & + write_restart_snow use ice_restart_driver, only: dumpfile use ice_restoring, only: restore_ice, ice_HaloRestore use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_horiz, step_dyn_ridge, step_radiation, & - biogeochemistry, save_init, step_dyn_wave + biogeochemistry, save_init, step_dyn_wave, step_snow use ice_timers, only: ice_timer_start, ice_timer_stop, & timer_diags, timer_column, timer_thermo, timer_bound, & timer_hist, timer_readwrite @@ -144,19 +145,28 @@ subroutine ice_step offset ! d(age)/dt time offset logical (kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_fsd, & + tr_iage, tr_FY, tr_lvl, tr_fsd, tr_snow, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo, tr_brine, tr_iso, tr_aero, & calc_Tsfc, skl_bgc, solve_zsal, z_tracers, wave_spec character(len=*), parameter :: subname = '(ice_step)' + character (len=char_len) :: plabeld + + if (debug_model) then + plabeld = 'beginning time step' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo + endif + call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc, skl_bgc_out=skl_bgc, & solve_zsal_out=solve_zsal, z_tracers_out=z_tracers, ktherm_out=ktherm, & wave_spec_out=wave_spec) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & tr_pond_topo_out=tr_pond_topo, tr_brine_out=tr_brine, tr_aero_out=tr_aero, & - tr_iso_out=tr_iso, tr_fsd_out=tr_fsd) + tr_iso_out=tr_iso, tr_fsd_out=tr_fsd, tr_snow_out=tr_snow) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -201,15 +211,33 @@ subroutine ice_step !----------------------------------------------------------------- if (calc_Tsfc) call prep_radiation (iblk) + if (debug_model) then + plabeld = 'post prep_radiation' + call debug_ice (iblk, plabeld) + endif !----------------------------------------------------------------- ! thermodynamics and biogeochemistry !----------------------------------------------------------------- call step_therm1 (dt, iblk) ! vertical thermodynamics + if (debug_model) then + plabeld = 'post step_therm1' + call debug_ice (iblk, plabeld) + endif + call biogeochemistry (dt, iblk) ! biogeochemistry + if (debug_model) then + plabeld = 'post biogeochemistry' + call debug_ice (iblk, plabeld) + endif + if (.not.prescribed_ice) & call step_therm2 (dt, iblk) ! ice thickness distribution thermo + if (debug_model) then + plabeld = 'post step_therm2' + call debug_ice (iblk, plabeld) + endif endif ! ktherm > 0 @@ -237,6 +265,12 @@ subroutine ice_step ! momentum, stress, transport call step_dyn_horiz (dt_dyn) + if (debug_model) then + plabeld = 'post step_dyn_horiz' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo ! iblk + endif ! ridging !$OMP PARALLEL DO PRIVATE(iblk) @@ -244,12 +278,24 @@ subroutine ice_step if (kridge > 0) call step_dyn_ridge (dt_dyn, ndtd, iblk) enddo !$OMP END PARALLEL DO + if (debug_model) then + plabeld = 'post step_dyn_ridge' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo ! iblk + endif ! clean up, update tendency diagnostics offset = c0 call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset) enddo + if (debug_model) then + plabeld = 'post dynamics' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo + endif endif ! not prescribed ice @@ -260,18 +306,36 @@ subroutine ice_step call ice_timer_start(timer_column) ! column physics call ice_timer_start(timer_thermo) ! thermodynamics + !----------------------------------------------------------------- + ! snow redistribution and metamorphosis + !----------------------------------------------------------------- + + if (tr_snow) then ! advanced snow physics + do iblk = 1, nblocks + call step_snow (dt, iblk) + enddo + call update_state (dt) ! clean up + endif + !MHRI: CHECK THIS OMP !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks if (ktherm >= 0) call step_radiation (dt, iblk) + if (debug_model) then + plabeld = 'post step_radiation' + call debug_ice (iblk, plabeld) + endif !----------------------------------------------------------------- ! get ready for coupling and the next time step !----------------------------------------------------------------- call coupling_prep (iblk) - + if (debug_model) then + plabeld = 'post coupling_prep' + call debug_ice (iblk, plabeld) + endif enddo ! iblk !$OMP END PARALLEL DO @@ -309,6 +373,7 @@ subroutine ice_step if (tr_pond_cesm) call write_restart_pond_cesm if (tr_pond_lvl) call write_restart_pond_lvl if (tr_pond_topo) call write_restart_pond_topo + if (tr_snow) call write_restart_snow if (tr_fsd) call write_restart_fsd if (tr_iso) call write_restart_iso if (tr_aero) call write_restart_aero @@ -634,11 +699,12 @@ subroutine sfcflux_to_ocn(nx_block, ny_block, & real (kind=dbl_kind) :: & puny, & ! + Lsub, & ! rLsub ! 1/Lsub character(len=*), parameter :: subname = '(sfcflux_to_ocn)' - call icepack_query_parameters(puny_out=puny) + call icepack_query_parameters(puny_out=puny, Lsub_out=Lsub) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) diff --git a/cicecore/drivers/nuopc/dmi/CICE_InitMod.F90 b/cicecore/drivers/nuopc/dmi/CICE_InitMod.F90 index 625348863..82f0ff0e8 100644 --- a/cicecore/drivers/nuopc/dmi/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/dmi/CICE_InitMod.F90 @@ -18,6 +18,7 @@ module CICE_InitMod use icepack_intfc, only: icepack_aggregate use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_init_snow use icepack_intfc, only: icepack_configure use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, & @@ -81,15 +82,14 @@ subroutine cice_init(mpi_comm) use ice_flux, only: init_coupler_flux, init_history_therm, & init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, & - get_forcing_atmo, get_forcing_ocn, get_wave_spec + get_forcing_atmo, get_forcing_ocn, get_wave_spec, init_snowtable use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & faero_default, faero_optics, alloc_forcing_bgc, fiso_default use ice_grid, only: init_grid1, init_grid2, alloc_grid use ice_history, only: init_hist, accum_hist use ice_restart_shared, only: restart, runtype use ice_init, only: input_data, init_state - use ice_init_column, only: init_thermo_vertical, init_shortwave, & - init_zbgc, input_zbgc, count_tracers + use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers use ice_kinds_mod use ice_restoring, only: ice_HaloRestore_init use ice_timers, only: timer_total, init_ice_timers, ice_timer_start @@ -99,7 +99,8 @@ subroutine cice_init(mpi_comm) mpi_comm ! communicator for sequential ccsm logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, & - tr_iso, tr_fsd, wave_spec + tr_iso, tr_fsd, wave_spec, tr_snow + character(len=char_len) :: snw_aging_table character(len=*), parameter :: subname = '(cice_init)' if (present(mpi_comm)) then @@ -176,7 +177,7 @@ subroutine cice_init(mpi_comm) call ice_HaloRestore_init ! restored boundary conditions call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & - wave_spec_out=wave_spec) + wave_spec_out=wave_spec, snw_aging_table_out=snw_aging_table) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(trim(subname), & file=__FILE__,line= __LINE__) @@ -190,7 +191,7 @@ subroutine cice_init(mpi_comm) call calc_timesteps ! update timestep counter if not using npt_unit="1" call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) - call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_query_tracer_flags(tr_iso_out=tr_iso, tr_snow_out=tr_snow) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(trim(subname), & file=__FILE__,line= __LINE__) @@ -222,8 +223,20 @@ subroutine cice_init(mpi_comm) call get_forcing_ocn(dt) ! ocean forcing from data #endif + ! snow aging lookup table initialization + if (tr_snow) then ! advanced snow physics + call icepack_init_snow() + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + if (snw_aging_table(1:4) /= 'test') then + call init_snowtable() + endif + endif + ! isotopes if (tr_iso) call fiso_default ! default values + ! aerosols ! if (tr_aero) call faero_data ! data file ! if (tr_zaero) call fzaero_data ! data file (gx1) @@ -252,12 +265,12 @@ subroutine init_restart use ice_calendar, only: calendar use ice_constants, only: c0 use ice_domain, only: nblocks - use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd, nslyr use ice_dyn_eap, only: read_restart_eap use ice_dyn_shared, only: kdyn use ice_grid, only: tmask use ice_init, only: ice_ic - use ice_init_column, only: init_age, init_FY, init_lvl, & + use ice_init_column, only: init_age, init_FY, init_lvl, init_snowtracers, & init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd use ice_restart_column, only: restart_age, read_restart_age, & @@ -265,6 +278,7 @@ subroutine init_restart restart_pond_cesm, read_restart_pond_cesm, & restart_pond_lvl, read_restart_pond_lvl, & restart_pond_topo, read_restart_pond_topo, & + restart_snow, read_restart_snow, & restart_fsd, read_restart_fsd, & restart_iso, read_restart_iso, & restart_aero, read_restart_aero, & @@ -279,12 +293,13 @@ subroutine init_restart iblk ! block index logical(kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & - tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + tr_pond_topo, tr_snow, tr_fsd, tr_iso, tr_aero, tr_brine, & skl_bgc, z_tracers, solve_zsal integer(kind=int_kind) :: & ntrcr integer(kind=int_kind) :: & nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_smice, nt_smliq, nt_rhos, nt_rsnw, & nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice character(len=*), parameter :: subname = '(init_restart)' @@ -299,10 +314,12 @@ subroutine init_restart call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & - tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + tr_snow_out=tr_snow, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & + nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw, & nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & @@ -399,6 +416,22 @@ subroutine init_restart enddo ! iblk endif ! .not. restart_pond endif + + ! snow redistribution/metamorphism + if (tr_snow) then + if (trim(runtype) == 'continue') restart_snow = .true. + if (restart_snow) then + call read_restart_snow + else + do iblk = 1, nblocks + call init_snowtracers(trcrn(:,:,nt_smice:nt_smice+nslyr-1,:,iblk), & + trcrn(:,:,nt_smliq:nt_smliq+nslyr-1,:,iblk), & + trcrn(:,:,nt_rhos :nt_rhos +nslyr-1,:,iblk), & + trcrn(:,:,nt_rsnw :nt_rsnw +nslyr-1,:,iblk)) + enddo ! iblk + endif + endif + ! floe size distribution if (tr_fsd) then if (trim(runtype) == 'continue') restart_fsd = .true. @@ -415,7 +448,7 @@ subroutine init_restart if (restart_iso) then call read_restart_iso else - do iblk = 1, nblocks + do iblk = 1, nblocks call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) enddo ! iblk diff --git a/cicecore/drivers/nuopc/dmi/CICE_RunMod.F90 b/cicecore/drivers/nuopc/dmi/CICE_RunMod.F90 index cfd519146..2d3e22973 100644 --- a/cicecore/drivers/nuopc/dmi/CICE_RunMod.F90 +++ b/cicecore/drivers/nuopc/dmi/CICE_RunMod.F90 @@ -41,7 +41,7 @@ module CICE_RunMod ! Philip W. Jones, LANL ! William H. Lipscomb, LANL - subroutine CICE_Run + subroutine CICE_Run(stop_now_cpl) use ice_calendar, only: istep, istep1, dt, stop_now, advance_timestep use ice_forcing, only: get_forcing_atmo, get_forcing_ocn, & @@ -54,6 +54,7 @@ subroutine CICE_Run logical (kind=log_kind) :: & tr_iso, tr_aero, tr_zaero, skl_bgc, z_tracers, wave_spec, tr_fsd character(len=*), parameter :: subname = '(CICE_Run)' + logical (kind=log_kind), optional, intent(in) :: stop_now_cpl !-------------------------------------------------------------------- ! initialize error code and step timer @@ -84,6 +85,9 @@ subroutine CICE_Run call advance_timestep() ! advance time + if (present(stop_now_cpl)) then + if (stop_now_cpl) return + endif #ifndef CICE_IN_NEMO #ifndef CICE_DMI if (stop_now >= 1) exit timeLoop @@ -142,7 +146,7 @@ subroutine ice_step use ice_boundary, only: ice_HaloUpdate use ice_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep - use ice_diagnostics, only: init_mass_diags, runtime_diags + use ice_diagnostics, only: init_mass_diags, runtime_diags, debug_model, debug_ice use ice_diagnostics_bgc, only: hbrine_diags, zsal_diags, bgc_diags use ice_domain, only: halo_info, nblocks use ice_dyn_eap, only: write_restart_eap @@ -155,12 +159,13 @@ subroutine ice_step use ice_restart_column, only: write_restart_age, write_restart_FY, & write_restart_lvl, write_restart_pond_cesm, write_restart_pond_lvl, & write_restart_pond_topo, write_restart_aero, write_restart_fsd, & - write_restart_iso, write_restart_bgc, write_restart_hbrine + write_restart_iso, write_restart_bgc, write_restart_hbrine, & + write_restart_snow use ice_restart_driver, only: dumpfile use ice_restoring, only: restore_ice, ice_HaloRestore use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_horiz, step_dyn_ridge, step_radiation, & - biogeochemistry, save_init, step_dyn_wave + biogeochemistry, save_init, step_dyn_wave, step_snow use ice_timers, only: ice_timer_start, ice_timer_stop, & timer_diags, timer_column, timer_thermo, timer_bound, & timer_hist, timer_readwrite @@ -174,19 +179,28 @@ subroutine ice_step offset ! d(age)/dt time offset logical (kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_fsd, & + tr_iage, tr_FY, tr_lvl, tr_fsd, tr_snow, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo, tr_brine, tr_iso, tr_aero, & calc_Tsfc, skl_bgc, solve_zsal, z_tracers, wave_spec character(len=*), parameter :: subname = '(ice_step)' + character (len=char_len) :: plabeld + + if (debug_model) then + plabeld = 'beginning time step' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo + endif + call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc, skl_bgc_out=skl_bgc, & solve_zsal_out=solve_zsal, z_tracers_out=z_tracers, ktherm_out=ktherm, & wave_spec_out=wave_spec) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & tr_pond_topo_out=tr_pond_topo, tr_brine_out=tr_brine, tr_aero_out=tr_aero, & - tr_iso_out=tr_iso, tr_fsd_out=tr_fsd) + tr_iso_out=tr_iso, tr_fsd_out=tr_fsd, tr_snow_out=tr_snow) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -223,14 +237,36 @@ subroutine ice_step if (calc_Tsfc) call prep_radiation (iblk) + if (debug_model) then + plabeld = 'post prep_radiation' + call debug_ice (iblk, plabeld) + endif + !----------------------------------------------------------------- ! thermodynamics and biogeochemistry !----------------------------------------------------------------- call step_therm1 (dt, iblk) ! vertical thermodynamics + + if (debug_model) then + plabeld = 'post step_therm1' + call debug_ice (iblk, plabeld) + endif + call biogeochemistry (dt, iblk) ! biogeochemistry + + if (debug_model) then + plabeld = 'post biogeochemistry' + call debug_ice (iblk, plabeld) + endif + call step_therm2 (dt, iblk) ! ice thickness distribution thermo + if (debug_model) then + plabeld = 'post step_therm2' + call debug_ice (iblk, plabeld) + endif + endif ! ktherm > 0 enddo ! iblk @@ -256,6 +292,13 @@ subroutine ice_step ! momentum, stress, transport call step_dyn_horiz (dt_dyn) + if (debug_model) then + plabeld = 'post step_dyn_horiz' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo ! iblk + endif + ! ridging !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks @@ -263,31 +306,66 @@ subroutine ice_step enddo !$OMP END PARALLEL DO + if (debug_model) then + plabeld = 'post step_dyn_ridge' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo ! iblk + endif + ! clean up, update tendency diagnostics offset = c0 call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset) enddo - !----------------------------------------------------------------- - ! albedo, shortwave radiation - !----------------------------------------------------------------- + if (debug_model) then + plabeld = 'post dynamics' + do iblk = 1, nblocks + call debug_ice (iblk, plabeld) + enddo + endif call ice_timer_start(timer_column) ! column physics call ice_timer_start(timer_thermo) ! thermodynamics + !----------------------------------------------------------------- + ! snow redistribution and metamorphosis + !----------------------------------------------------------------- + + if (tr_snow) then ! advanced snow physics + do iblk = 1, nblocks + call step_snow (dt, iblk) + enddo + call update_state (dt) ! clean up + endif + !MHRI: CHECK THIS OMP !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks + !----------------------------------------------------------------- + ! albedo, shortwave radiation + !----------------------------------------------------------------- + if (ktherm >= 0) call step_radiation (dt, iblk) + if (debug_model) then + plabeld = 'post step_radiation' + call debug_ice (iblk, plabeld) + endif + !----------------------------------------------------------------- ! get ready for coupling and the next time step !----------------------------------------------------------------- call coupling_prep (iblk) + if (debug_model) then + plabeld = 'post coupling_prep' + call debug_ice (iblk, plabeld) + endif + enddo ! iblk !$OMP END PARALLEL DO @@ -325,6 +403,7 @@ subroutine ice_step if (tr_pond_cesm) call write_restart_pond_cesm if (tr_pond_lvl) call write_restart_pond_lvl if (tr_pond_topo) call write_restart_pond_topo + if (tr_snow) call write_restart_snow if (tr_fsd) call write_restart_fsd if (tr_iso) call write_restart_iso if (tr_aero) call write_restart_aero diff --git a/configuration/scripts/Makefile b/configuration/scripts/Makefile index 51c36cee3..07ee380a8 100644 --- a/configuration/scripts/Makefile +++ b/configuration/scripts/Makefile @@ -73,7 +73,6 @@ RM := rm AR := ar .SUFFIXES: -.SUFFIXES: .F90 .F .c .o .PHONY: all cice libcice targets target db_files db_flags clean realclean helloworld calchk sumchk bcstchk all: $(EXEC) @@ -167,13 +166,13 @@ libcice: $(OBJS) @ echo "$(AR) -r $(EXEC) $(OBJS)" $(AR) -r $(EXEC) $(OBJS) -.c.o: +%.o : %.c $(CC) $(CFLAGS) $(CPPDEFS) $(INCLDIR) $< -.F.o: +%.o : %.F %.d $(FC) -c $(FFLAGS) $(FIXEDFLAGS) $(CPPDEFS) $(INCLDIR) $< -.F90.o: +%.o : %.F90 %.d $(FC) -c $(FFLAGS) $(FREEFLAGS) $(CPPDEFS) $(MODDIR) $(INCLDIR) $< clean: diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 3dec72963..bb44663eb 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -141,7 +141,8 @@ Cstar = 20 Cf = 17. Ktens = 0. - e_ratio = 2. + e_yieldcurve = 2. + e_plasticpot = 2. seabed_stress = .false. seabed_stress_method = 'LKD' k1 = 7.5 @@ -223,7 +224,7 @@ &forcing_nml formdrag = .false. - atmbndy = 'default' + atmbndy = 'similarity' rotate_wind = .true. calc_strair = .true. calc_Tsfc = .true. diff --git a/configuration/scripts/machines/Macros.daley_intel b/configuration/scripts/machines/Macros.daley_intel index a434ffdb3..3755946f7 100644 --- a/configuration/scripts/machines/Macros.daley_intel +++ b/configuration/scripts/machines/Macros.daley_intel @@ -13,8 +13,9 @@ FFLAGS := -fp-model source -convert big_endian -assume byterecl -ftz -traceb #-xHost ifeq ($(ICE_BLDDEBUG), true) - FFLAGS += -O0 -g -check -fpe0 -ftrapuv -fp-model except -check noarg_temp_created -init=snan,arrays -# -heap-arrays 1024 + FFLAGS += -O0 -g -check -fpe0 -ftrapuv -fp-model except -check noarg_temp_created +#-heap-arrays 1024 +#-init=snan,arrays else FFLAGS += -O2 endif diff --git a/configuration/scripts/machines/env.cheyenne_gnu b/configuration/scripts/machines/env.cheyenne_gnu index 3bfe59c31..8ddc443b1 100755 --- a/configuration/scripts/machines/env.cheyenne_gnu +++ b/configuration/scripts/machines/env.cheyenne_gnu @@ -29,6 +29,13 @@ if ($ICE_IOTYPE =~ pio*) then endif endif +if ($?ICE_TEST) then +if ($ICE_TEST =~ qcchk*) then + module load python + source /glade/u/apps/opt/ncar_pylib/ncar_pylib.csh default +endif +endif + endif limit coredumpsize unlimited diff --git a/configuration/scripts/machines/env.cheyenne_intel b/configuration/scripts/machines/env.cheyenne_intel index 4a430622e..28df6647d 100755 --- a/configuration/scripts/machines/env.cheyenne_intel +++ b/configuration/scripts/machines/env.cheyenne_intel @@ -29,6 +29,13 @@ if ($ICE_IOTYPE =~ pio*) then endif endif +if ($?ICE_TEST) then +if ($ICE_TEST =~ qcchk*) then + module load python + source /glade/u/apps/opt/ncar_pylib/ncar_pylib.csh default +endif +endif + endif limit coredumpsize unlimited diff --git a/configuration/scripts/machines/env.cheyenne_pgi b/configuration/scripts/machines/env.cheyenne_pgi index 693692842..d492129fb 100755 --- a/configuration/scripts/machines/env.cheyenne_pgi +++ b/configuration/scripts/machines/env.cheyenne_pgi @@ -29,6 +29,13 @@ if ($ICE_IOTYPE =~ pio*) then endif endif +if ($?ICE_TEST) then +if ($ICE_TEST =~ qcchk*) then + module load python + source /glade/u/apps/opt/ncar_pylib/ncar_pylib.csh default +endif +endif + endif limit coredumpsize unlimited diff --git a/configuration/scripts/options/set_env.box2001 b/configuration/scripts/options/set_env.box2001 deleted file mode 100644 index a3f7c10f5..000000000 --- a/configuration/scripts/options/set_env.box2001 +++ /dev/null @@ -1 +0,0 @@ -setenv NICELYR 1 diff --git a/configuration/scripts/options/set_nml.alt03 b/configuration/scripts/options/set_nml.alt03 index 255d77261..98e794735 100644 --- a/configuration/scripts/options/set_nml.alt03 +++ b/configuration/scripts/options/set_nml.alt03 @@ -19,7 +19,7 @@ sw_dtemp = 0.02d0 tfrz_option = 'linear_salt' revised_evp = .false. Ktens = 0. -e_ratio = 2. +e_yieldcurve = 2. seabed_stress = .true. use_bathymetry = .true. l_mpond_fresh = .true. diff --git a/configuration/scripts/options/set_nml.atmbndyconstant b/configuration/scripts/options/set_nml.atmbndyconstant new file mode 100644 index 000000000..8e2a68192 --- /dev/null +++ b/configuration/scripts/options/set_nml.atmbndyconstant @@ -0,0 +1 @@ +atmbndy = 'constant' diff --git a/configuration/scripts/options/set_nml.atmbndymixed b/configuration/scripts/options/set_nml.atmbndymixed new file mode 100644 index 000000000..c79a26fb5 --- /dev/null +++ b/configuration/scripts/options/set_nml.atmbndymixed @@ -0,0 +1 @@ +atmbndy = 'mixed' diff --git a/configuration/scripts/options/set_nml.gx1coreii b/configuration/scripts/options/set_nml.gx1coreii deleted file mode 100644 index 13b8db59e..000000000 --- a/configuration/scripts/options/set_nml.gx1coreii +++ /dev/null @@ -1,10 +0,0 @@ -year_init = 1997 -use_leap_years = .false. -use_restart_time = .true. -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' -fyear_init = 2005 -ycycle = 1 -atm_data_format = 'bin' -atm_data_type = 'LYq' -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/COREII' -precip_units = 'mm_per_sec' diff --git a/configuration/scripts/options/set_nml.qc_nonbfb b/configuration/scripts/options/set_nml.qc_nonbfb deleted file mode 100644 index e34f75c35..000000000 --- a/configuration/scripts/options/set_nml.qc_nonbfb +++ /dev/null @@ -1,8 +0,0 @@ -npt = 87600 -dt = 1800.0 -dumpfreq = 'm' -dumpfreq_n = 12 -diagfreq = 24 -histfreq = 'd','x','x','x','x' -f_hi = 'd' -hist_avg = .false. diff --git a/configuration/scripts/options/set_nml.qcnonbfb b/configuration/scripts/options/set_nml.qcnonbfb new file mode 100644 index 000000000..a965b863c --- /dev/null +++ b/configuration/scripts/options/set_nml.qcnonbfb @@ -0,0 +1,16 @@ +dt = 3456.0 +npt_unit = 'y' +npt = 5 +year_init = 2005 +month_init = 1 +day_init = 1 +sec_init = 0 +use_leap_years = .false. +fyear_init = 2005 +ycycle = 1 +dumpfreq = 'm' +dumpfreq_n = 12 +diagfreq = 24 +histfreq = 'd','x','x','x','x' +f_hi = 'd' +hist_avg = .false. diff --git a/configuration/scripts/options/set_nml.snwITDrdg b/configuration/scripts/options/set_nml.snwitdrdg similarity index 100% rename from configuration/scripts/options/set_nml.snwITDrdg rename to configuration/scripts/options/set_nml.snwitdrdg diff --git a/configuration/scripts/parse_namelist.sh b/configuration/scripts/parse_namelist.sh index ea539a2d0..dcb0d1ccc 100755 --- a/configuration/scripts/parse_namelist.sh +++ b/configuration/scripts/parse_namelist.sh @@ -10,7 +10,7 @@ filename=$1 filemods=$2 #echo "$0 $1 $2" -echo "running parse_namelist.sh" +echo "running ${scriptname}" foundstring="FoundSTRING" vnamearray=() valuearray=() @@ -43,11 +43,9 @@ do fi done - #sed -i 's|\(^\s*'"$vname"'\s*\=\s*\)\(.*$\)|\1'"$value"'|g' $filename - cp ${filename} ${filename}.check - sed -i.sedbak -e 's|\(^[[:space:]]*'"$vname"'[[:space:]]*=[[:space:]]*\)\(.*$\)|\1'"$foundstring"'|g' ${filename}.check - grep -q ${foundstring} ${filename}.check - if [ $? -eq 0 ]; then + grep -q "^[[:space:]]*${vname}[[:space:]]*=" $filename + grepout=$? + if [ ${grepout} -eq 0 ]; then sed -i.sedbak -e 's|\(^[[:space:]]*'"$vname"'[[:space:]]*=[[:space:]]*\)\(.*$\)|\1'"$value"'|g' ${filename} if [[ "${found}" == "${foundstring}" ]]; then vnamearray+=($vname) @@ -55,17 +53,17 @@ do else valuearray[$found]=${value} fi - if [[ -e "${filename}.sedbak" ]]; then - rm ${filename}.sedbak - fi else echo "${scriptname} ERROR: parsing error for ${vname}" exit -99 fi - rm ${filename}.check ${filename}.check.sedbak fi done < "$filemods" +if [[ -e "${filename}.sedbak" ]]; then + rm ${filename}.sedbak +fi + exit 0 diff --git a/configuration/scripts/parse_namelist_from_env.sh b/configuration/scripts/parse_namelist_from_env.sh index 4d829450f..4c25d358d 100755 --- a/configuration/scripts/parse_namelist_from_env.sh +++ b/configuration/scripts/parse_namelist_from_env.sh @@ -5,10 +5,11 @@ if [[ "$#" -ne 1 ]]; then exit -1 fi +scriptname=`basename "$0"` filename=$1 #echo "$0 $1" -echo "running parse_namelist_from_env.sh" +echo "running $scriptname" sed -i.sedbak -e 's|ICE_SANDBOX|'"${ICE_SANDBOX}"'|g' $filename sed -i.sedbak -e 's|ICE_MACHINE_INPUTDATA|'"${ICE_MACHINE_INPUTDATA}"'|g' $filename diff --git a/configuration/scripts/parse_settings.sh b/configuration/scripts/parse_settings.sh index d6ed31c15..a3f432801 100755 --- a/configuration/scripts/parse_settings.sh +++ b/configuration/scripts/parse_settings.sh @@ -10,7 +10,7 @@ filename=$1 filemods=$2 #echo "$0 $1 $2" -echo "running parse_settings.sh" +echo "running ${scriptname}" foundstring="FoundSTRING" vnamearray=() valuearray=() @@ -23,8 +23,11 @@ do else #vname=`echo $line | sed "s|\(^\s*set\S*\)\s\{1,100\}\(\S*\)\s\{1,100\}\(\S*\).*$|\2|g"` #value=`echo $line | sed "s|\(^\s*set\S*\)\s\{1,100\}\(\S*\)\s\{1,100\}\(\S*\).*$|\3|g"` - vname=`echo $line | sed "s|\(^[[:space:]]*set[^[:space:]]*\)[[:space:]][[:space:]]*\([^[:space:]]*\)[[:space:]][[:space:]]*\([^[:space:]]*\).*$|\2|g"` + vname=`echo $line | sed "s|\(^[[:space:]]*set[^[:space:]]*\)[[:space:]][[:space:]]*\([^[:space:]]*\).*$|\2|g"` value=`echo $line | sed "s|\(^[[:space:]]*set[^[:space:]]*\)[[:space:]][[:space:]]*\([^[:space:]]*\)[[:space:]][[:space:]]*\([^[:space:]]*\).*$|\3|g"` + if [[ "${value}" == "${line}" ]]; then + value="" + fi # echo "$line $vname $value" found=${foundstring} @@ -43,22 +46,27 @@ do fi done - #sed -i 's|\(^\s*set.* '"$vname"' \)[^#]*\(#*.*$\)|\1 '"$value"' \2|g' $filename - sed -i.sedbak -e 's|\(^[[:space:]]*set.* '"$vname"' \)[^#]*\(#*.*$\)|\1 '"$value"' \2|g' $filename - - if [[ "${found}" == "${foundstring}" ]]; then - vnamearray+=($vname) - valuearray+=($value) + grep -q "^[[:space:]]*set.* ${vname}[[:space:]]*" $filename + grepout=$? + if [ ${grepout} -eq 0 ]; then + sed -i.sedbak -e 's|\(^[[:space:]]*set.* '"$vname"' \)[^#]*\(#*.*$\)|\1 '"$value"' \2|g' $filename + if [[ "${found}" == "${foundstring}" ]]; then + vnamearray+=($vname) + valuearray+=($value) + else + valuearray[$found]=${value} + fi else - valuearray[$found]=${value} - fi - - if [[ -e "${filename}.sedbak" ]]; then - rm ${filename}.sedbak + echo "${scriptname} ERROR: parsing error for ${vname}" + exit -99 fi fi done < "$filemods" +if [[ -e "${filename}.sedbak" ]]; then + rm ${filename}.sedbak +fi + exit 0 diff --git a/configuration/scripts/tests/QC/cice.t-test.py b/configuration/scripts/tests/QC/cice.t-test.py index 6f2c7e89b..c84583baa 100755 --- a/configuration/scripts/tests/QC/cice.t-test.py +++ b/configuration/scripts/tests/QC/cice.t-test.py @@ -177,7 +177,10 @@ def stage_one(data_d, num_files, mean_d, variance_d): df = n_eff - 1 # Read in t_crit table - nfid = nc.Dataset("configuration/scripts/tests/QC/CICE_t_critical_p0.8.nc", 'r') + if os.path.exists('./CICE_t_critical_p0.8.nc'): + nfid = nc.Dataset("./CICE_t_critical_p0.8.nc", 'r') + else: + nfid = nc.Dataset("configuration/scripts/tests/QC/CICE_t_critical_p0.8.nc", 'r') df_table = nfid.variables['df'][:] t_crit_table = nfid.variables['tcrit'][:] nfid.close() @@ -238,7 +241,10 @@ def stage_one(data_d, num_files, mean_d, variance_d): t_val = mean_d / np.sqrt(variance_d / num_files) # Find t_crit from the nearest value on the Lookup Table Test - nfid = nc.Dataset("configuration/scripts/tests/QC/CICE_Lookup_Table_p0.8_n1825.nc", 'r') + if os.path.exists('./CICE_Lookup_Table_p0.8_n1825.nc'): + nfid = nc.Dataset("./CICE_Lookup_Table_p0.8_n1825.nc", 'r') + else: + nfid = nc.Dataset("configuration/scripts/tests/QC/CICE_Lookup_Table_p0.8_n1825.nc", 'r') r1_table = nfid.variables['r1'][:] t_crit_table = nfid.variables['tcrit'][:] nfid.close() diff --git a/configuration/scripts/tests/QC/gen_qc_cases.csh b/configuration/scripts/tests/QC/gen_qc_cases.csh index 7a395f1c2..7d2db73b9 100755 --- a/configuration/scripts/tests/QC/gen_qc_cases.csh +++ b/configuration/scripts/tests/QC/gen_qc_cases.csh @@ -183,9 +183,9 @@ endif # Generate the non-BFB but non-climate-changing case echo "Generating nonbfb case" if ($testid != $spval) then - set result = `./cice.setup $options -s qc_nonbfb,long --testid qc_test_$testid | grep 'Test case dir\|already exists'` + set result = `./cice.setup $options -s qcnonbfb,long --testid qc_test_$testid | grep 'Test case dir\|already exists'` else - set result = `./cice.setup $options -s qc_nonbfb,long --testid qc_test | grep 'Test case dir\|already exists'` + set result = `./cice.setup $options -s qcnonbfb,long --testid qc_test | grep 'Test case dir\|already exists'` endif set nonbfb_dir = `echo "$result" | awk '{print$NF}'` if ($nonbfb_dir == "exists") then diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 4da4dd110..e4c376ad4 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -58,9 +58,9 @@ restart gx3 4x2 fsd12,debug,short smoke gx3 8x2 fsd12ww3,diag24,run1day smoke gx3 4x1 isotope,debug restart gx3 8x2 isotope -smoke gx3 4x1 snwITDrdg,snwgrain,icdefault,debug +smoke gx3 4x1 snwitdrdg,snwgrain,icdefault,debug smoke gx3 4x1 snw30percent,icdefault,debug -restart gx3 8x2 snwITDrdg,icdefault,snwgrain +restart gx3 8x2 snwitdrdg,icdefault,snwgrain restart gx3 4x4 gx3ncarbulk,iobinary restart gx3 4x4 histall,precision8,cdf64 smoke gx3 30x1 bgcz,histall @@ -70,5 +70,6 @@ smoke gx3 8x2 diag24,run5day,zsal,debug restart gx3 8x2 zsal restart gx3 8x2 gx3ncarbulk,debug restart gx3 4x4 gx3ncarbulk,diag1 -restart gx1 24x1 gx1coreii,short smoke gx3 4x1 calcdragio +restart gx3 4x2 atmbndyconstant +restart gx3 4x2 atmbndymixed diff --git a/configuration/scripts/tests/baseline.script b/configuration/scripts/tests/baseline.script index 6f13807e3..ac69d49a0 100644 --- a/configuration/scripts/tests/baseline.script +++ b/configuration/scripts/tests/baseline.script @@ -36,6 +36,12 @@ if (${ICE_BASECOM} != ${ICE_SPVAL}) then ${ICE_CASEDIR}/casescripts/comparelog.csh ${base_file} ${test_file} notcicefile set bfbstatus = $status + else if (${ICE_TEST} =~ qcchk*) then + set test_dir = ${ICE_RUNDIR} + set base_dir = ${ICE_BASELINE}/${ICE_BASECOM}/${ICE_TESTNAME} + ${ICE_SANDBOX}/configuration/scripts/tests/QC/cice.t-test.py ${base_dir} ${test_dir} + set bfbstatus = $status + else set test_dir = ${ICE_RUNDIR}/restart @@ -116,6 +122,35 @@ endif if (${ICE_BFBCOMP} != ${ICE_SPVAL}) then + echo "PEND ${ICE_TESTNAME} bfbcomp ${ICE_BFBCOMP}" >> ${ICE_CASEDIR}/test_output + if (${ICE_BFBCOMP} != ${ICE_TESTNAME}) then + # Check if the baseline job is complete + set job = `grep " ${ICE_BFBCOMP} " ../suite.jobs | sed 's|^[^0-9]*\([0-9]*\).*$|\1|g'` + echo "checking on Job $job" + set qstatjob = 1 + set cnt = 0 + if (${job} =~ [0-9]*) then + while ($qstatjob) + ${ICE_MACHINE_QSTAT} $job >&/dev/null + set qstatus = $status +# echo $job $qstatus + if ($qstatus != 0) then + echo "Job $job completed" + set qstatjob = 0 + else + @ cnt = $cnt + 1 + echo "Waiting for $job to complete $cnt" + sleep 60 # Sleep for 1 minute, so as not to overwhelm the queue manager + if ($cnt > 100) then + echo "No longer waiting for $job to complete" + set qstatjob = 0 # Abandon check after 100 sleep 60 checks + endif + endif +# echo $qstatjob + end + endif + endif + if (${ICE_TEST} == "logbfb") then set test_file = `ls -1t ${ICE_RUNDIR}/cice.runlog* | head -1` set base_file = `ls -1t ${ICE_RUNDIR}/../${ICE_BFBCOMP}.${ICE_TESTID}/cice.runlog* | head -1` @@ -127,6 +162,16 @@ if (${ICE_BFBCOMP} != ${ICE_SPVAL}) then ${ICE_CASEDIR}/casescripts/comparelog.csh ${base_file} ${test_file} set bfbstatus = $status + + else if (${ICE_TEST} =~ qcchk*) then + set test_dir = ${ICE_RUNDIR} + set base_dir = ${ICE_RUNDIR}/../${ICE_BFBCOMP}.${ICE_TESTID} + ${ICE_SANDBOX}/configuration/scripts/tests/QC/cice.t-test.py ${base_dir} ${test_dir} + set bfbstatus = $status + # expecting failure, so switch value + if (${ICE_TEST} =~ qcchkf*) then + @ bfbstatus = 1 - $bfbstatus + endif else set test_dir = ${ICE_RUNDIR}/restart set base_dir = ${ICE_RUNDIR}/../${ICE_BFBCOMP}.${ICE_TESTID}/restart @@ -140,6 +185,9 @@ if (${ICE_BFBCOMP} != ${ICE_SPVAL}) then set bfbstatus = $status endif + mv -f ${ICE_CASEDIR}/test_output ${ICE_CASEDIR}/test_output.prev + cat ${ICE_CASEDIR}/test_output.prev | grep -iv "${ICE_TESTNAME} bfbcomp " >! ${ICE_CASEDIR}/test_output + rm -f ${ICE_CASEDIR}/test_output.prev if (${bfbstatus} == 0) then echo "PASS ${ICE_TESTNAME} bfbcomp ${ICE_BFBCOMP}" >> ${ICE_CASEDIR}/test_output echo "bfb baseline and test dataset are identical" diff --git a/configuration/scripts/tests/create_fails.csh b/configuration/scripts/tests/create_fails.csh new file mode 100755 index 000000000..d22fe01ad --- /dev/null +++ b/configuration/scripts/tests/create_fails.csh @@ -0,0 +1,24 @@ +#!/bin/csh + +echo " " +set tmpfile = create_fails.tmp +set outfile = fails.ts + +./results.csh >& /dev/null +cat results.log | grep ' run\| test' | grep -v "#" | grep -v PASS | cut -f 2 -d " " | sort -u >! $tmpfile + +echo "# Test Grid PEs Sets" >! $outfile +foreach line ( "`cat $tmpfile`" ) + #echo $line + set test = `echo $line | cut -d "_" -f 3` + set grid = `echo $line | cut -d "_" -f 4` + set pes = `echo $line | cut -d "_" -f 5` + set opts = `echo $line | cut -d "_" -f 6- | sed 's/_/,/g'` + echo "$test $grid $pes $opts" >> $outfile +end + +rm $tmpfile +echo "$0 done" +echo "Failed tests can be rerun with the test suite file...... $outfile" +echo "To run a new test suite, copy $outfile to the top directory and do something like" +echo " ./cice.setup --suite $outfile ..." diff --git a/configuration/scripts/tests/decomp_suite.ts b/configuration/scripts/tests/decomp_suite.ts index 9c82c5d27..c39c3ddfe 100644 --- a/configuration/scripts/tests/decomp_suite.ts +++ b/configuration/scripts/tests/decomp_suite.ts @@ -3,7 +3,6 @@ restart gx3 4x2x25x29x4 dslenderX2 restart gx1 64x1x16x16x10 dwghtfile restart gbox180 16x1x6x6x60 dspacecurve,debugblocks decomp gx3 4x2x25x29x5 none -sleep 30 restart gx3 1x1x50x58x4 droundrobin,thread restart_gx3_4x2x25x29x4_dslenderX2 restart gx3 4x1x25x116x1 dslenderX1,thread restart_gx3_4x2x25x29x4_dslenderX2 restart gx3 6x2x4x29x18 dspacecurve restart_gx3_4x2x25x29x4_dslenderX2 @@ -27,7 +26,6 @@ restart gx3 8x1x25x29x4 drakeX2,thread restart_gx3_4x2x25x2 smoke gx3 4x2x25x29x4 debug,run2day,dslenderX2 smoke gx1 64x1x16x16x10 debug,run2day,dwghtfile smoke gbox180 16x1x6x6x60 debug,run2day,dspacecurve,debugblocks -sleep 30 smoke gx3 1x1x25x58x8 debug,run2day,droundrobin,thread smoke_gx3_4x2x25x29x4_debug_dslenderX2_run2day smoke gx3 20x1x5x116x1 debug,run2day,dslenderX1,thread smoke_gx3_4x2x25x29x4_debug_dslenderX2_run2day smoke gx3 6x2x4x29x18 debug,run2day,dspacecurve smoke_gx3_4x2x25x29x4_debug_dslenderX2_run2day diff --git a/configuration/scripts/tests/first_suite.ts b/configuration/scripts/tests/first_suite.ts index bf6c813f6..31eba9fb7 100644 --- a/configuration/scripts/tests/first_suite.ts +++ b/configuration/scripts/tests/first_suite.ts @@ -1,5 +1,6 @@ # Test Grid PEs Sets BFB-compare smoke gx3 8x2 diag1,run5day restart gx3 4x2x25x29x4 dslenderX2 +smoke gx3 4x2x25x29x4 debug,run2day,dslenderX2 logbfb gx3 4x2x25x29x4 dslenderX2,diag1,reprosum smoke gx3 1x2 run2day diff --git a/configuration/scripts/tests/nothread_suite.ts b/configuration/scripts/tests/nothread_suite.ts index da1267e86..616741aa2 100644 --- a/configuration/scripts/tests/nothread_suite.ts +++ b/configuration/scripts/tests/nothread_suite.ts @@ -19,7 +19,7 @@ restart gx3 16x1 gx3ncarbulk,iobinary restart gx3 12x1 alt01 restart gx3 16x1 alt02 restart gx3 8x1 alt03 -restart gx3 16x1 alt04 +restart gx3 16x1x5x29x6 alt04 restart gx3 16x1 alt05 restart gx3 20x1 alt06 restart gx3 18x1 alt01,debug,short @@ -66,7 +66,7 @@ restart gx3 16x1x8x10x10 droundrobin restart_gx3_8x1x25x29x2_ restart gx3 6x1x50x58x1 droundrobin restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 8x1x19x19x5 droundrobin restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 20x1x5x29x20 dsectrobin,short restart_gx3_8x1x25x29x2_dslenderX2 -restart gx3 32x1x5x10x10 drakeX2 restart_gx3_8x1x25x29x2_dslenderX2 +restart gx3 32x1x5x10x12 drakeX2 restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 16x1x8x10x10 droundrobin,maskhalo restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 4x1x25x29x4 droundrobin restart_gx3_8x1x25x29x2_dslenderX2 diff --git a/configuration/scripts/tests/prod_suite.ts b/configuration/scripts/tests/prod_suite.ts index 8793dfed2..04982adb1 100644 --- a/configuration/scripts/tests/prod_suite.ts +++ b/configuration/scripts/tests/prod_suite.ts @@ -1,4 +1,6 @@ # Test Grid PEs Sets BFB-compare -smoke gx1 64x1 qc,medium -smoke gx1 64x2 gx1prod,long,run10year - +qcchk gx3 72x1 qc,medium qcchk_gx3_72x1_medium_qc +qcchk gx1 144x1 qc,medium +smoke gx1 144x2 gx1prod,long,run10year +qcchkf gx3 72x1 qc,medium,alt02 qcchk_gx3_72x1_medium_qc +qcchk gx3 72x1 qcnonbfb,medium qcchk_gx3_72x1_medium_qc diff --git a/configuration/scripts/tests/reprosum_suite.ts b/configuration/scripts/tests/reprosum_suite.ts index dd6a6d56b..a7f3fe5bc 100644 --- a/configuration/scripts/tests/reprosum_suite.ts +++ b/configuration/scripts/tests/reprosum_suite.ts @@ -1,7 +1,6 @@ # Test Grid PEs Sets BFB-compare logbfb gx3 4x2x25x29x4 dslenderX2,diag1,reprosum #logbfb gx3 4x2x25x29x4 dslenderX2,diag1 -sleep 60 logbfb gx3 1x1x50x58x4 droundrobin,diag1,thread,maskhalo,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum logbfb gx3 4x1x25x116x1 dslenderX1,diag1,thread,maskhalo,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum logbfb gx3 1x20x5x29x80 dsectrobin,diag1,short,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum diff --git a/configuration/scripts/tests/test_qcchk.script b/configuration/scripts/tests/test_qcchk.script new file mode 100644 index 000000000..81b5f05fc --- /dev/null +++ b/configuration/scripts/tests/test_qcchk.script @@ -0,0 +1,36 @@ + +cp ${ICE_SANDBOX}/configuration/scripts/tests/QC/CICE_t_critical_p0.8.nc . +cp ${ICE_SANDBOX}/configuration/scripts/tests/QC/CICE_Lookup_Table_p0.8_n1825.nc . + +#---------------------------------------------------- +# Run the CICE model +# cice.run returns -1 if run did not complete successfully + +./cice.run +set res="$status" + +set log_file = `ls -t1 ${ICE_RUNDIR}/cice.runlog* | head -1` +set ttimeloop = `grep TimeLoop ${log_file} | grep Timer | cut -c 22-32` +set tdynamics = `grep Dynamics ${log_file} | grep Timer | cut -c 22-32` +set tcolumn = `grep Column ${log_file} | grep Timer | cut -c 22-32` +if (${ttimeloop} == "") set ttimeloop = -1 +if (${tdynamics} == "") set tdynamics = -1 +if (${tcolumn} == "") set tcolumn = -1 + +mv -f ${ICE_CASEDIR}/test_output ${ICE_CASEDIR}/test_output.prev +cat ${ICE_CASEDIR}/test_output.prev | grep -iv "${ICE_TESTNAME} run" >! ${ICE_CASEDIR}/test_output +mv -f ${ICE_CASEDIR}/test_output ${ICE_CASEDIR}/test_output.prev +cat ${ICE_CASEDIR}/test_output.prev | grep -iv "${ICE_TESTNAME} test" >! ${ICE_CASEDIR}/test_output +rm -f ${ICE_CASEDIR}/test_output.prev + +set grade = PASS +if ( $res != 0 ) then + set grade = FAIL + echo "$grade ${ICE_TESTNAME} run ${ttimeloop} ${tdynamics} ${tcolumn}" >> ${ICE_CASEDIR}/test_output + echo "$grade ${ICE_TESTNAME} test " >> ${ICE_CASEDIR}/test_output + exit 99 +endif + +echo "$grade ${ICE_TESTNAME} run ${ttimeloop} ${tdynamics} ${tcolumn}" >> ${ICE_CASEDIR}/test_output +echo "$grade ${ICE_TESTNAME} test " >> ${ICE_CASEDIR}/test_output + diff --git a/configuration/scripts/tests/test_qcchkf.script b/configuration/scripts/tests/test_qcchkf.script new file mode 100644 index 000000000..81b5f05fc --- /dev/null +++ b/configuration/scripts/tests/test_qcchkf.script @@ -0,0 +1,36 @@ + +cp ${ICE_SANDBOX}/configuration/scripts/tests/QC/CICE_t_critical_p0.8.nc . +cp ${ICE_SANDBOX}/configuration/scripts/tests/QC/CICE_Lookup_Table_p0.8_n1825.nc . + +#---------------------------------------------------- +# Run the CICE model +# cice.run returns -1 if run did not complete successfully + +./cice.run +set res="$status" + +set log_file = `ls -t1 ${ICE_RUNDIR}/cice.runlog* | head -1` +set ttimeloop = `grep TimeLoop ${log_file} | grep Timer | cut -c 22-32` +set tdynamics = `grep Dynamics ${log_file} | grep Timer | cut -c 22-32` +set tcolumn = `grep Column ${log_file} | grep Timer | cut -c 22-32` +if (${ttimeloop} == "") set ttimeloop = -1 +if (${tdynamics} == "") set tdynamics = -1 +if (${tcolumn} == "") set tcolumn = -1 + +mv -f ${ICE_CASEDIR}/test_output ${ICE_CASEDIR}/test_output.prev +cat ${ICE_CASEDIR}/test_output.prev | grep -iv "${ICE_TESTNAME} run" >! ${ICE_CASEDIR}/test_output +mv -f ${ICE_CASEDIR}/test_output ${ICE_CASEDIR}/test_output.prev +cat ${ICE_CASEDIR}/test_output.prev | grep -iv "${ICE_TESTNAME} test" >! ${ICE_CASEDIR}/test_output +rm -f ${ICE_CASEDIR}/test_output.prev + +set grade = PASS +if ( $res != 0 ) then + set grade = FAIL + echo "$grade ${ICE_TESTNAME} run ${ttimeloop} ${tdynamics} ${tcolumn}" >> ${ICE_CASEDIR}/test_output + echo "$grade ${ICE_TESTNAME} test " >> ${ICE_CASEDIR}/test_output + exit 99 +endif + +echo "$grade ${ICE_TESTNAME} run ${ttimeloop} ${tdynamics} ${tcolumn}" >> ${ICE_CASEDIR}/test_output +echo "$grade ${ICE_TESTNAME} test " >> ${ICE_CASEDIR}/test_output + diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 0a04b5e26..85acbece3 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -65,7 +65,7 @@ either Celsius or Kelvin units). "atm_data_dir", "directory for atmospheric forcing data", "" "atm_data_format", "format of atmospheric forcing files", "" "atm_data_type", "type of atmospheric forcing", "" - "atmbndy", "atmo boundary layer parameterization (‘default’ or ‘constant’)", "" + "atmbndy", "atmo boundary layer parameterization ('similarity', ‘constant’, or 'mixed')", "" "avail_hist_fields", "type for history field data", "" "awtidf", "weighting factor for near-ir, diffuse albedo", "0.36218" "awtidr", "weighting factor for near-ir, direct albedo", "0.00182" @@ -202,9 +202,12 @@ either Celsius or Kelvin units). "eps13", "a small number", "10\ :math:`^{-13}`" "eps16", "a small number", "10\ :math:`^{-16}`" "esno(n)", "energy of melting of snow per unit area (in category n)", "J/m\ :math:`^2`" + "etax2", "2 x eta (shear viscous coefficient)", "kg/s" "evap", "evaporative water flux", "kg/m\ :math:`^2`/s" "ew_boundary_type", "type of east-west boundary condition", "" - "eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36" + "eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36" + "e_yieldcurve", "yield curve minor/major axis ratio", "2" + "e_plasticpot", "plastic potential minor/major axis ratio", "2" "**F**", "", "" "faero_atm", "aerosol deposition rate", "kg/m\ :math:`^2`/s" "faero_ocn", "aerosol flux to the ocean", "kg/m\ :math:`^2`/s" @@ -515,7 +518,6 @@ either Celsius or Kelvin units). "print_global", "if true, print global data", "F" "print_points", "if true, print point data", "F" "processor_shape", "descriptor for processor aspect ratio", "" - "prs_sig", "replacement pressure", "N/m" "Pstar", "ice strength parameter", "2.75\ :math:`\times`\ 10\ :math:`^4`\ N/m\ :math:`^2`" "puny", "a small positive number", "1\ :math:`\times`\ 10\ :math:`^{-11}`" "**Q**", "", "" @@ -540,6 +542,7 @@ either Celsius or Kelvin units). "rdg_shear", "shear for ridging", "1/s" "real_kind", "definition of single precision real", "selected_real_kind(6)" "refindx", "refractive index of sea ice", "1.310" + "rep_prs", "replacement pressure", "N/m" "revp", "real(revised_evp)", "" "restart", "if true, initialize ice state from file", "T" "restart_age", "if true, read age restart file", "" @@ -734,6 +737,7 @@ either Celsius or Kelvin units). "yieldstress11(12, 22)", "yield stress tensor components", "" "year_init", "the initial year", "" "**Z**", "", "" + "zetax2", "2 x zeta (bulk viscous coefficient)", "kg/s" "zlvl", "atmospheric level height (momentum)", "m" "zlvs", "atmospheric level height (scalars)", "m" "zref", "reference height for stability", "10. m" diff --git a/doc/source/developer_guide/dg_forcing.rst b/doc/source/developer_guide/dg_forcing.rst index aea6d8ef6..d3c406b43 100644 --- a/doc/source/developer_guide/dg_forcing.rst +++ b/doc/source/developer_guide/dg_forcing.rst @@ -150,21 +150,6 @@ Users are encouraged to switch to the JRA55 (see :ref:`JRA55forcing`) dataset. atmosphere forcing dataset may be deprecated in the future. -.. _LYqforcing: - -LYq Atmosphere Forcing -------------------------- - -The LYq (:cite:`Hunke07`) forcing was used in earlier standalone -runs on the gx1 grid, and the -Consortium continues to do some very limited testing with this forcing dataset. -This dataset is largely based on the CORE II data. -Monthly average data for cldf and fsnow is read while 6-hourly data for Qa, Tair, -uatm, and vatm are read with other fields derived or set by default. -Users are encouraged to switch to the JRA55 (see :ref:`JRA55forcing`) dataset. This -atmosphere forcing dataset may be deprecated in the future. - - .. _defaultforcing: Default Atmosphere Forcing diff --git a/doc/source/developer_guide/dg_scripts.rst b/doc/source/developer_guide/dg_scripts.rst index 50853b3ea..dac5e9a52 100644 --- a/doc/source/developer_guide/dg_scripts.rst +++ b/doc/source/developer_guide/dg_scripts.rst @@ -78,7 +78,8 @@ are the three scripts that modify **ice_in** and **cice.settings**. To add new options, just add new files to the **configurations/scripts/options/** directory with appropriate names and syntax. The set_nml file syntax is the same as namelist syntax and the set_env files are consistent with csh setenv syntax. See other files for -examples of the syntax. +examples of the syntax. The name of the option (i.e. diag1, debug, bgcISPOL) should not +have any special characters in the name as this can impact scripts usage. .. _build: @@ -163,7 +164,7 @@ it's working properly. .. _dev_validation: -Code Validation Script +QC Process Validation ---------------------- The code validation (aka QC or quality control) test validates non bit-for-bit model changes. The directory @@ -193,9 +194,9 @@ to the ``cice.setup`` script. These options include: * ``--queue`` : Queue for the batch submission * ``--testid`` : test ID, user-defined id for testing -The script creates 4 test cases, with testIDs ``qc_base``, ``qc_bfb``, ``qc_nonbfb``, +The script creates 4 test cases, with testIDs ``qc_base``, ``qc_bfb``, ``qc_test``, and ``qc_fail``. ``qc_base`` is the base test case with the default QC namelist. -``qc_bfb`` is identical to ``qc_base``. ``qc_nonbfb`` is a test that is not bit-for-bit +``qc_bfb`` is identical to ``qc_base``. ``qc_test`` is a test that is not bit-for-bit when compared to ``qc_base``, but not climate changing. ``qc_fail`` is a test that is not bit-for-bit and also climate changing. @@ -222,7 +223,7 @@ To perform the QC validation, execute the following commands. # From the CICE base directory cp configuration/scripts/tests/QC/gen_qc_cases.csh . - cp configuration/scripts/tests/QC/compare_qc_cases.csh + cp configuration/scripts/tests/QC/compare_qc_cases.csh . # Create the required test cases ./gen_qc_cases.csh -m --acct diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 7c2a45a35..e4afa0c46 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -60,7 +60,7 @@ @string{CRST @string{IJHPCA={Int. J High Perform. Comput. Appl}} @string{PTRSA={Philos. Trans. Royal Soc. A}} @string{SIAMJCP={SIAM J. Sci. Comput.}} - +@string{TC={The Cryosphere}} % ********************************************** @@ -979,6 +979,17 @@ @article{Roach18 volume = {123}, year = {2018} } + +@Article{Ringeisen21 + author = "D. Ringeisen and L.B. Tremblay and M. Losch", + title = "{Non-normal flow rules affect fracture angles in sea ice viscous-plastic rheologies}", + journal = TC, + year = {2021}, + volume = {15}, + pages = {2873-2888}, + url = {https://doi.org/10.5194/tc-15-2873-2021} +} + @Article{Tsujino18, author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.", title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}", diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index d4e209d8a..9c529b8ec 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -350,9 +350,9 @@ Following again the LKD method, the seabed stress coefficients are finally expre .. _internal-stress: -*************** +******************************** Internal stress -*************** +******************************** For convenience we formulate the stress tensor :math:`\bf \sigma` in terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, @@ -378,15 +378,6 @@ CICE can output the internal ice pressure which is an important field to support The internal ice pressure (``sigP``) is the average of the normal stresses multiplied by :math:`-1` and is therefore simply equal to :math:`-\sigma_1/2`. -Following the approach of :cite:`Konig10` (see also :cite:`Lemieux16`), the -elliptical yield curve can be modified such that the ice has isotropic tensile strength. -The tensile strength :math:`T_p` is expressed as a fraction of the ice strength :math:`P`, that is :math:`T_p=k_t P` -where :math:`k_t` should be set to a value between 0 and 1 (this can -be changed at runtime with the namelist parameter ``Ktens``). The ice -strength :math:`P` is a function of the ice thickness distribution as -described in the `Icepack -Documentation`_. - .. _stress-vp: Viscous-Plastic @@ -395,28 +386,45 @@ Viscous-Plastic The VP constitutive law is given by .. math:: - \sigma_{ij} = 2 \eta \dot{\epsilon}_{ij} + (\zeta - \eta) D_D - P_R(1 - k_t)\frac{\delta_{ij}}{2} + \sigma_{ij} = 2 \eta \dot{\epsilon}_{ij} + (\zeta - \eta) D_D - P_R\frac{\delta_{ij}}{2} :label: vp-const -where :math:`\eta` and :math:`\zeta` are the bulk and shear viscosities. +where :math:`\eta` and :math:`\zeta` are the bulk and shear viscosities and +:math:`P_R` is a “replacement pressure” (see :cite:`Geiger98`, for example), +which serves to prevent residual ice motion due to spatial +variations of the ice strength :math:`P` when the strain rates are exactly zero. + An elliptical yield curve is used, with the viscosities given by .. math:: \zeta = {P(1+k_t)\over 2\Delta}, .. math:: - \eta = {P(1+k_t)\over {2\Delta e^2}}, + \eta = e_g^{-2} \zeta, where .. math:: - \Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2} + \Delta = \left[D_D^2 + {e_f^2\over e_g^4}\left(D_T^2 + D_S^2\right)\right]^{1/2}. + +The ice strength :math:`P` is a function of the ice thickness distribution as +described in the `Icepack Documentation `_. + +Two modifications to the standard VP rheology of :cite:`Hibler79` are available. +First, following the approach of :cite:`Konig10` (see also :cite:`Lemieux16`), the +elliptical yield curve can be modified such that the ice has isotropic tensile strength. +The tensile strength is expressed as a fraction of :math:`P`, that is :math:`k_t P` +where :math:`k_t` should be set to a value between 0 and 1 (this can +be changed at runtime with the namelist parameter ``Ktens``). -and :math:`P_R` is a “replacement pressure” (see :cite:`Geiger98`, for -example), which serves to prevent residual ice motion due to spatial -variations of :math:`P` when the rates of strain are exactly zero. +Second, while :math:`e_f` is the ratio of the major and minor axes of the elliptical yield curve, the parameter +:math:`e_g` characterizes the plastic potential, i.e. another ellipse that decouples the flow rule from the +yield curve (:cite:`Ringeisen21`). :math:`e_f` and :math:`e_g` are respectively called ``e_yieldcurve`` and ``e_plasticpot`` in the code and +can be set in the namelist. The plastic potential can lead to more realistic fracture angles between linear kinematic features. :cite:`Ringeisen21` suggest to set :math:`e_f` to a value larger than 1 and to have :math:`e_g < e_f`. -The parameter :math:`e` is the ratio of the major and minor axes of the elliptical yield curve, also called the ellipse aspect ratio. It can be changed using the namelist parameter ``e_ratio``. +By default, the namelist parameters are set to :math:`e_f=e_g=2` and :math:`k_t=0` which correspond to the standard VP rheology. + +There are four options in the code for solving the sea ice momentum equation with a VP formulation: the standard EVP approach, a 1d EVP solver, the revised EVP approach and an implicit Picard solver. The modifications to the yield curve and to the flow rule described above are available for these four different solution methods. .. _stress-evp: @@ -428,7 +436,7 @@ regularized version of the VP constitutive law :eq:`vp-const`. The constitutive .. math:: {1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} - + {P_R(1-k_t)\over 2\zeta} = D_D, \\ + + {P_R\over 2\zeta} = D_D, \\ :label: sig1 .. math:: @@ -455,32 +463,103 @@ parameter less than one. Including the modification proposed by :cite:`Bouillon1 .. math:: \begin{aligned} {\partial\sigma_1\over\partial t} + {\sigma_1\over 2T} - + {P_R(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta} D_D, \\ - {\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {P(1+k_t)\over - 2Te^2\Delta} D_T,\\ + + {P_R\over 2T} &=& {\zeta \over T} D_D, \\ + {\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {\eta \over + T} D_T,\\ {\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 2T} &=& - {P(1+k_t)\over 4Te^2\Delta}D_S.\end{aligned} + {\eta \over 2T}D_S.\end{aligned} Once discretized in time, these last three equations are written as .. math:: \begin{aligned} {(\sigma_1^{k+1}-\sigma_1^{k})\over\Delta t_e} + {\sigma_1^{k+1}\over 2T} - + {P_R^k(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta^k} D_D^k, \\ - {(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {P(1+k_t)\over - 2Te^2\Delta^k} D_T^k,\\ + + {P_R^k\over 2T} &=& {\zeta^k\over T} D_D^k, \\ + {(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {\eta^k \over + T} D_T^k,\\ {(\sigma_{12}^{k+1}-\sigma_{12}^{k})\over\Delta t_e} + {\sigma_{12}^{k+1}\over 2T} &=& - {P(1+k_t)\over 4Te^2\Delta^k}D_S^k,\end{aligned} + {\eta^k \over 2T}D_S^k,\end{aligned} :label: sigdisc where :math:`k` denotes again the subcycling step. All coefficients on the left-hand side are constant except for :math:`P_R`. This modification compensates for the decreased efficiency of including -the viscosity terms in the subcycling. (Note that the viscosities do not -appear explicitly.) Choices of the parameters used to define :math:`E`, +the viscosity terms in the subcycling. Choices of the parameters used to define :math:`E`, :math:`T` and :math:`\Delta t_e` are discussed in Sections :ref:`revp` and :ref:`parameters`. +.. _evp1d: + +1d EVP solver +~~~~~~~~~~~~~ + +The standard EVP solver iterates hundreds of times, where each iteration includes a communication through MPI and a limited number of calculations. This limits how much the solver can be optimized as the speed is primarily determined by the communication. The 1d EVP solver avoids the communication by utilizing shared memory, which removes the requirement for calls to the MPI communicator. As a consequence of this the potential scalability of the code is improved. The performance is best on shared memory but the solver is also functional on MPI and hybrid MPI/OpenMP setups as it will run on the master processor alone. + +The scalability of geophysical models is in general terms limited by the memory usage. In order to optimize this the 1d EVP solver solves the same equations that are outlined in the section :ref:`stress-evp` but it transforms all matrices to vectors (1d matrices) as this compiles better with the computer hardware. The vectorization and the contiguous placement of arrays in the memory makes it easier for the compiler to optimize the code and pass pointers instead of copying the vectors. The 1d solver is not supported for tripole grids and the code will abort if this combination is attempted. + +.. _revp: + +Revised EVP approach +~~~~~~~~~~~~~~~~~~~~ + +The revised EVP approach is based on a pseudo-time iterative scheme :cite:`Lemieux12`, :cite:`Bouillon13`, :cite:`Kimmritz15`. By construction, the revised EVP approach should lead to the VP solution +(given the right numerical parameters and a sufficiently large number of iterations). To do so, the inertial term is formulated such that it matches the backward Euler approach of +implicit solvers and there is an additional term for the pseudo-time iteration. Hence, with the revised approach, the discretized momentum equations :eq:`umom` and :eq:`vmom` become + +.. math:: + {\beta^*(u^{k+1}-u^k)\over\Delta t_e} + {m(u^{k+1}-u^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)} u^{k+1} + - {\left(mf+{\tt vrel}\sin\theta\right)} v^{k+1} + = {{\partial\sigma_{1j}^{k+1}\over\partial x_j}} + + {\tau_{ax} - mg{\partial H_\circ\over\partial x} } + + {\tt vrel} {\left(U_w\cos\theta-V_w\sin\theta\right)}, + :label: umomr + +.. math:: + {\beta^*(v^{k+1}-v^k)\over\Delta t_e} + {m(v^{k+1}-v^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)}v^{k+1} + + {\left(mf+{\tt vrel}\sin\theta\right)} u^{k+1} + = {{\partial\sigma_{2j}^{k+1}\over\partial x_j}} + + {\tau_{ay} - mg{\partial H_\circ\over\partial y} } + + {\tt vrel}{\left(U_w\sin\theta+V_w\cos\theta\right)}, + :label: vmomr + +where :math:`\beta^*` is a numerical parameter and :math:`u^n, v^n` are the components of the previous time level solution. +With :math:`\beta=\beta^* \Delta t \left( m \Delta t_e \right)^{-1}` :cite:`Bouillon13`, these equations can be written as + +.. math:: + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta\ + C_b \right)}_{\tt cca} u^{k+1} + - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} + = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} + + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} + + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t}(\beta u^k + u^n), + :label: umomr2 + +.. math:: + \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} + + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta + C_b \right)}_{\tt cca}v^{k+1} + = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} + + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} + + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t}(\beta v^k + v^n), + :label: vmomr2 + +At this point, the solutions :math:`u^{k+1}` and :math:`v^{k+1}` are obtained in the same manner as for the standard EVP approach (see equations :eq:`cevpuhat` to :eq:`cevpb`). + +Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite:`Bouillon13`, the stress equations in :eq:`sigdisc` become + +.. math:: + \begin{aligned} + {\alpha (\sigma_1^{k+1}-\sigma_1^{k})} + {\sigma_1^{k}} + + {P_R^k} &=& 2 \zeta^k D_D^k, \\ + {\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& 2 \eta^k D_T^k,\\ + {\alpha (\sigma_{12}^{k+1}-\sigma_{12}^{k})} + {\sigma_{12}^{k}} &=& + \eta^k D_S^k,\end{aligned} + +where as opposed to the classic EVP, the second term in each equation is at iteration :math:`k` :cite:`Bouillon13`. Also, as opposed to the classic EVP, +:math:`\Delta t_e` times the number of subcycles (or iterations) does not need to be equal to the advective time step :math:`\Delta t`. +Finally, as with the classic EVP approach, the stresses are initialized using the previous time level values. +The revised EVP is activated by setting the namelist parameter ``revised_evp = true``. +In the code :math:`\alpha` is ``arlx`` and :math:`\beta` is ``brlx``. The values of ``arlx`` and ``brlx`` can be set in the namelist. +It is recommended to use large values of these parameters and to set :math:`\alpha=\beta` :cite:`Kimmritz15`. + .. _stress-eap: Elastic-Anisotropic-Plastic @@ -668,78 +747,3 @@ rheology we compute the area loss rate due to ridging as Both ridging rate and sea ice strength are computed in the outer loop of the dynamics. - -.. _revp: - -**************** -Revised approach -**************** - -The revised EVP approach is based on a pseudo-time iterative scheme :cite:`Lemieux12`, :cite:`Bouillon13`, :cite:`Kimmritz15`. By construction, the revised EVP approach should lead to the VP solution -(given the right numerical parameters and a sufficiently large number of iterations). To do so, the inertial term is formulated such that it matches the backward Euler approach of -implicit solvers and there is an additional term for the pseudo-time iteration. Hence, with the revised approach, the discretized momentum equations :eq:`umom` and :eq:`vmom` become - -.. math:: - {\beta^*(u^{k+1}-u^k)\over\Delta t_e} + {m(u^{k+1}-u^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)} u^{k+1} - - {\left(mf+{\tt vrel}\sin\theta\right)} v^{k+1} - = {{\partial\sigma_{1j}^{k+1}\over\partial x_j}} - + {\tau_{ax} - mg{\partial H_\circ\over\partial x} } - + {\tt vrel} {\left(U_w\cos\theta-V_w\sin\theta\right)}, - :label: umomr - -.. math:: - {\beta^*(v^{k+1}-v^k)\over\Delta t_e} + {m(v^{k+1}-v^n)\over\Delta t} + {\left({\tt vrel} \cos\theta + C_b \right)}v^{k+1} - + {\left(mf+{\tt vrel}\sin\theta\right)} u^{k+1} - = {{\partial\sigma_{2j}^{k+1}\over\partial x_j}} - + {\tau_{ay} - mg{\partial H_\circ\over\partial y} } - + {\tt vrel}{\left(U_w\sin\theta+V_w\cos\theta\right)}, - :label: vmomr - -where :math:`\beta^*` is a numerical parameter and :math:`u^n, v^n` are the components of the previous time level solution. -With :math:`\beta=\beta^* \Delta t \left( m \Delta t_e \right)^{-1}` :cite:`Bouillon13`, these equations can be written as - -.. math:: - \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta\ + C_b \right)}_{\tt cca} u^{k+1} - - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} - = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} - + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} - + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t}(\beta u^k + u^n), - :label: umomr2 - -.. math:: - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} - + \underbrace{\left((\beta+1){m\over\Delta t}+{\tt vrel} \cos\theta + C_b \right)}_{\tt cca}v^{k+1} - = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} - + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} - + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t}(\beta v^k + v^n), - :label: vmomr2 - -At this point, the solutions :math:`u^{k+1}` and :math:`v^{k+1}` are obtained in the same manner as for the standard EVP approach (see equations :eq:`cevpuhat` to :eq:`cevpb`). - -Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite:`Bouillon13`, the stress equations in :eq:`sigdisc` become - -.. math:: - \begin{aligned} - {\alpha (\sigma_1^{k+1}-\sigma_1^{k})} + {\sigma_1^{k}} - + {P_R^k(1-k_t)} &=& {P(1+k_t)\over \Delta^k} D_D^k, \\ - {\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& {P(1+k_t)\over - e^2\Delta^k} D_T^k,\\ - {\alpha (\sigma_{12}^{k+1}-\sigma_{12}^{k})} + {\sigma_{12}^{k}} &=& - {P(1+k_t)\over 2e^2\Delta^k}D_S^k,\end{aligned} - -where as opposed to the classic EVP, the second term in each equation is at iteration :math:`k` :cite:`Bouillon13`. Also, as opposed to the classic EVP, -:math:`\Delta t_e` times the number of subcycles (or iterations) does not need to be equal to the advective time step :math:`\Delta t`. -Finally, as with the classic EVP approach, the stresses are initialized using the previous time level values. -The revised EVP is activated by setting the namelist parameter ``revised_evp = true``. -In the code :math:`\alpha` is ``arlx`` and :math:`\beta` is ``brlx``. The values of ``arlx`` and ``brlx`` can be set in the namelist. -It is recommended to use large values of these parameters and to set :math:`\alpha=\beta` :cite:`Kimmritz15`. - -.. _evp1d: - -**************** -1d EVP solver -**************** - -The standard EVP solver iterates hundreds of times, where each iteration includes a communication through MPI and a limited number of calculations. This limits how much the solver can be optimized as the speed is primarily determined by the communication. The 1d EVP solver avoids the communication by utilizing shared memory, which removes the requirement for calls to the MPI communicator. As a consequence of this the potential scalability of the code is improved. The performance is best on shared memory but the solver is also functional on MPI and hybrid MPI/OpenMP setups as it will run on the master processor alone. - -The scalability of geophysical models is in general terms limited by the memory usage. In order to optimize this the 1d EVP solver solves the same equations that are outlined in the section :ref:`stress-evp` but it transforms all matrices to vectors (1d matrices) as this compiles better with the computer hardware. The vectorization and the contiguous placement of arrays in the memory makes it easier for the compiler to optimize the code and pass pointers instead of copying the vectors. The 1d solver is not supported for tripole grids and the code will abort if this combination is attempted. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 65a50120b..23e6951fc 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -528,8 +528,10 @@ forcing_nml :widths: 15, 15, 30, 15 "", "", "", "" - "``atmbndy``", "``constant``", "bulk transfer coefficients", "``default``" - "", "``default``", "stability-based boundary layer", "" + "``atmbndy``", "string", "bulk transfer coefficients", "``similarity``" + "", "``similarity``", "stability-based boundary layer", "" + "", "``constant``", "constant-based boundary layer", "" + "", "``mixed``", "stability-based boundary layer for wind stress, constant-based for sensible+latent heat fluxes", "" "``atmiter_conv``", "real", "convergence criteria for ustar", "0.0" "``atm_data_dir``", "string", "path to atmospheric forcing data directory", "" "``atm_data_format``", "``bin``", "read direct access binary atmo forcing file format", "``bin``" @@ -540,7 +542,6 @@ forcing_nml "", "``JRA55_gx1``", "JRA55 forcing data for gx1 grid :cite:`Tsujino18`", "" "", "``JRA55_gx3``", "JRA55 forcing data for gx3 grid :cite:`Tsujino18`", "" "", "``JRA55_tx1``", "JRA55 forcing data for tx1 grid :cite:`Tsujino18`", "" - "", "``LYq``", "COREII Large-Yeager (AOMIP) forcing data :cite:`Large09`", "" "", "``monthly``", "monthly forcing data", "" "", "``ncar``", "NCAR bulk forcing data", "" "", "``oned``", "column forcing data", "" diff --git a/doc/source/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index ce5c2ef41..005f4f851 100644 --- a/doc/source/user_guide/ug_testing.rst +++ b/doc/source/user_guide/ug_testing.rst @@ -297,6 +297,13 @@ results.csh script in the testsuite.[testid]:: cd testsuite.[testid] ./results.csh +The script **create_fails.csh** will process the output from results.csh and generate a new +test suite file, **fails.ts**, from the failed tests. +**fails.ts** can then be edited and passed into ``cice.setup --suite fails.ts ...`` to rerun +subsets of failed tests to more efficiently move thru the development, testing, and +validation process. However, a full test suite should be run on the final development +version of the code. + To report the test results, as is required for Pull Requests to be accepted into the master the CICE Consortium code see :ref:`testreporting`. @@ -411,8 +418,10 @@ The *cice.setup** options ``--setup-only``, ``--setup-build``, and ``--setup-bui which means by default the test suite builds and submits the jobs. By defining other values for those environment variables, users can control the suite script. When using **suite.submit** manually, the string ``true`` (all lowercase) is the only string that will turn on a feature, and both SUITE_RUN and SUITE_SUBMIT cannot be true at the same time. -By leveraging the **cice.setup** command line arguments ``--setup-only``, ``--setup-build``, and ``--setup-build-run`` as well as the environment variables SUITE_BUILD, SUITE_RUN, and SUITE_SUBMIT, users can run **cice.setup** and **suite.submit** in various combinations to quickly setup, setup and build, submit, resubmit, run interactively, or rebuild and resubmit full testsuites quickly and easily. See below for an example. +By leveraging the **cice.setup** command line arguments ``--setup-only``, ``--setup-build``, and ``--setup-build-run`` as well as the environment variables SUITE_BUILD, SUITE_RUN, and SUITE_SUBMIT, users can run **cice.setup** and **suite.submit** in various combinations to quickly setup, setup and build, submit, resubmit, run interactively, or rebuild and resubmit full testsuites quickly and easily. See :ref:`examplesuites` for an example. + +.. _examplesuites: Test Suite Examples ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1121,7 +1130,7 @@ If the regression comparisons fail, then you may want to run the QC test, # From the updated sandbox # Generate the same test case(s) as the baseline using options or namelist changes to activate new code modifications - ./cice.setup -m onyx -e intel --test smoke -g gx1 -p 44x1 -testid qc_test -s qc,medium + ./cice.setup -m onyx -e intel --test smoke -g gx1 -p 44x1 --testid qc_test -s qc,medium cd onyx_intel_smoke_gx1_44x1_medium_qc.qc_test # modify ice_in to activate the namelist options that were determined above ./cice.build @@ -1130,7 +1139,8 @@ If the regression comparisons fail, then you may want to run the QC test, # Wait for runs to finish # Perform the QC test - cp configuration/scripts/tests/QC/cice.t-test.py + # From the updated sandbox + cp configuration/scripts/tests/QC/cice.t-test.py . ./cice.t-test.py /p/work/turner/CICE_RUNS/onyx_intel_smoke_gx1_44x1_medium_qc.qc_base \ /p/work/turner/CICE_RUNS/onyx_intel_smoke_gx1_44x1_medium_qc.qc_test diff --git a/doc/source/user_guide/ug_troubleshooting.rst b/doc/source/user_guide/ug_troubleshooting.rst index cd8f1acaf..d20b14ffc 100644 --- a/doc/source/user_guide/ug_troubleshooting.rst +++ b/doc/source/user_guide/ug_troubleshooting.rst @@ -77,7 +77,7 @@ using `runtype` = ‘initial’. Binary restart files that were provided with CICE v4.1 were made using the BL99 thermodynamics with 4 layers and 5 thickness categories (`kcatbound` = 0) and therefore can not be used for the default CICE v5 and later configuration (7 layers). In addition, CICE’s -default restart file format is now  instead of binary. +default restart file format is now NetCDF instead of binary. Restarting a run using `runtype` = ‘continue’ requires restart data for all tracers used in the new run. If tracer restart data is not @@ -126,7 +126,7 @@ conflicts in module dependencies. *print\_state* (**ice\_diagnostics.F90**) Print the ice state and forcing fields for a given grid cell. -`forcing\_diag` = true (**ice\_in**) +`debug\_forcing` = true (**ice\_in**) Print numerous diagnostic quantities associated with input forcing. `debug\_blocks` = true (**ice\_in**) diff --git a/icepack b/icepack index 34c8e688b..152bd701e 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 34c8e688bf7f3008cf84093cd703cf8cfe068eda +Subproject commit 152bd701e0cf3ec4385e5ce81918ba94e7a791cb