From 554a30278662cf68dc1e14c0197fdf66db88458c Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Fri, 10 Sep 2021 12:29:17 -0600 Subject: [PATCH 01/19] Add allocation checks for accumulating history variables. (#631) --- cicecore/cicedynB/analysis/ice_history.F90 | 58 +++++++++++++------ .../cicedynB/analysis/ice_history_bgc.F90 | 58 +++++++++---------- .../cicedynB/analysis/ice_history_drag.F90 | 3 + 3 files changed, 72 insertions(+), 47 deletions(-) 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 From 8c63ea0facbc584b0ec81e04df1dc25cff81aa4a Mon Sep 17 00:00:00 2001 From: TRasmussen <33480590+TillRasmussen@users.noreply.github.com> Date: Wed, 15 Sep 2021 09:58:31 +0200 Subject: [PATCH 02/19] initialized all eap history variables to zero on land (#632) --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 9c52bb888..9ecc79305 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 From 6c4160e77eebc5344a92d36fbcfe15e4ceef7b71 Mon Sep 17 00:00:00 2001 From: Mads Hvid Ribergaard <38077893+mhrib@users.noreply.github.com> Date: Mon, 27 Sep 2021 23:20:35 +0200 Subject: [PATCH 03/19] sensible+latent heatfluxes using linear bulk formula (#633) * 'heatflux_linear' flag: sensible+latent heatfluxes using traditional linear bulk formula * Add option atmbndy='mixed' boundary layer condition * New options for 'atmbndy' * For backward compability, rename "atmbndy" = default to similarity * Change wording for text related to atmbndy * Change wording for text related to atmbndy * Change wording for text related to atmbndy * Spelling error etc. * Update Icepack with required atmbndy changes * Added test for atmbndy={constant,mixed} * Add test for atmbndy={constant,mixed} * Test for atmbndy={constant,mixed} * clean up. Renamed without underscores --- cicecore/cicedynB/general/ice_init.F90 | 23 +++++++++++++++---- configuration/scripts/ice_in | 2 +- .../scripts/options/set_nml.atmbndyconstant | 1 + .../scripts/options/set_nml.atmbndymixed | 1 + configuration/scripts/tests/base_suite.ts | 2 ++ doc/source/cice_index.rst | 2 +- doc/source/user_guide/ug_case_settings.rst | 6 +++-- icepack | 2 +- 8 files changed, 29 insertions(+), 10 deletions(-) create mode 100644 configuration/scripts/options/set_nml.atmbndyconstant create mode 100644 configuration/scripts/options/set_nml.atmbndymixed diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 3d102217a..2e67af51c 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -448,7 +448,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 @@ -1214,6 +1214,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' @@ -1641,13 +1649,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/configuration/scripts/ice_in b/configuration/scripts/ice_in index 3dec72963..443ff1cbb 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -223,7 +223,7 @@ &forcing_nml formdrag = .false. - atmbndy = 'default' + atmbndy = 'similarity' rotate_wind = .true. calc_strair = .true. calc_Tsfc = .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/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 4da4dd110..2c75bffe1 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -72,3 +72,5 @@ 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/doc/source/cice_index.rst b/doc/source/cice_index.rst index 0a04b5e26..19428c87e 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" diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 65a50120b..2b4e26c72 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``" diff --git a/icepack b/icepack index 34c8e688b..f9c9e480f 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 34c8e688bf7f3008cf84093cd703cf8cfe068eda +Subproject commit f9c9e480f6ce482317734be80719178c8e1b5121 From dc18f61f9dbe36a9b3969e8e6f3ca47e287ace49 Mon Sep 17 00:00:00 2001 From: TRasmussen <33480590+TillRasmussen@users.noreply.github.com> Date: Wed, 6 Oct 2021 19:20:29 +0200 Subject: [PATCH 04/19] Reduce number of calls to calc_ffrac (#638) * initialized all eap history variables to zero on land * Step_a calls calc_fracv 6 times. This can be reduced to 3 as they 2 and two calcilates the same with limited difference * removed timers used for testing --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 94 +++++++++------------- 1 file changed, 38 insertions(+), 56 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 9ecc79305..83374d4dd 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -1939,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 @@ -2010,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 @@ -2042,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)) @@ -2068,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 From db96c724b921314afcd809283b16f295399ab682 Mon Sep 17 00:00:00 2001 From: Mads Hvid Ribergaard <38077893+mhrib@users.noreply.github.com> Date: Thu, 7 Oct 2021 23:24:31 +0200 Subject: [PATCH 05/19] Update for new snow scheme + introduce stop_now_cpl for coupling (#641) --- cicecore/drivers/nuopc/dmi/CICE_InitMod.F90 | 55 +++++++++--- cicecore/drivers/nuopc/dmi/CICE_RunMod.F90 | 97 +++++++++++++++++++-- 2 files changed, 132 insertions(+), 20 deletions(-) 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 From d794491747ecbe5df37540142d8e09b0dafe1e33 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Thu, 7 Oct 2021 14:30:42 -0700 Subject: [PATCH 06/19] Update README.md Change master to main --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 457714b8a..13a011213 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) +[![Travis-CI](https://travis-ci.org/CICE-Consortium/CICE.svg?branch=main)](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) From 9eddfc5aaa86536cb263230c6c3d29592e0ba9ac Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Thu, 7 Oct 2021 14:33:57 -0700 Subject: [PATCH 07/19] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 13a011213..45e4a6a1a 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![Travis-CI](https://travis-ci.org/CICE-Consortium/CICE.svg?branch=main)](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=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) From 9e4d244d7feb6ad9d03cc488c34c1f22c4ea7612 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Fri, 8 Oct 2021 18:24:05 -0400 Subject: [PATCH 08/19] Express rheology term as a function of viscous coefficients for EVP and rEVP (#639) * added subroutine to calc zeta and eta * new subroutine now also calculates the replacement pressure * stress calc in evp are now with zeta, eta and rep_prs * corrected minor compiling issues * improved comments and added references for papers * Modif to doc for new code with viscous coefficients * Minor modifs to doc * minor correction to QC testing doc * Small corrections following Philippes comments * Changed eta,zeta to etax2,zetax2 --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 100 ++++++++++-------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 72 ++++++++++++- .../scripts/machines/Macros.daley_intel | 5 +- doc/source/cice_index.rst | 4 +- doc/source/science_guide/sg_dynamics.rst | 30 +++--- doc/source/user_guide/ug_testing.rst | 5 +- 6 files changed, 149 insertions(+), 67 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 276c8bb79..775f2caf1 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 , & @@ -669,6 +672,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 +689,52 @@ 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, & + Deltase, Deltasw, & + zetax2ne, zetax2nw, & + zetax2se, zetax2sw, & + etax2ne, etax2nw, & + etax2se, etax2sw, & + rep_prsne, rep_prsnw, & + rep_prsse, rep_prssw ) + !----------------------------------------------------------------- ! 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..167bc6a2c 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 @@ -864,7 +865,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) ! @@ -1358,6 +1359,75 @@ subroutine strain_rates (nx_block, ny_block, & 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, & + Deltase, Deltasw, & + zetax2ne, zetax2nw, & + zetax2se, zetax2sw, & + etax2ne, etax2nw, & + etax2se, etax2sw, & + rep_prsne, rep_prsnw,& + rep_prsse, rep_prssw ) + + real (kind=dbl_kind), intent(in):: & + strength, tinyarea ! at the t-point + + real (kind=dbl_kind), intent(in):: & + Deltane, Deltanw, Deltase, Deltasw ! Delta at each corner + + real (kind=dbl_kind), intent(out):: & + zetax2ne, zetax2nw, zetax2se, zetax2sw, & ! zetax2 at each corner + etax2ne, etax2nw, etax2se, etax2sw, & ! etax2 at each corner + rep_prsne, rep_prsnw, rep_prsse, rep_prssw ! replacement pressure + + ! local variables + real (kind=dbl_kind) :: & + tmpcalc + + ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code + +! if (trim(yield_curve) == 'ellipse') then + + tmpcalc = strength/max(Deltane,tinyarea) ! northeast + zetax2ne = (c1+Ktens)*tmpcalc + rep_prsne = (c1-Ktens)*tmpcalc*Deltane + etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot + + tmpcalc = strength/max(Deltanw,tinyarea) ! northwest + zetax2nw = (c1+Ktens)*tmpcalc + rep_prsnw = (c1-Ktens)*tmpcalc*Deltanw + etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot + + tmpcalc = strength/max(Deltase,tinyarea) ! southeast + zetax2se = (c1+Ktens)*tmpcalc + rep_prsse = (c1-Ktens)*tmpcalc*Deltase + etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot + + tmpcalc = strength/max(Deltasw,tinyarea) ! southwest + zetax2sw = (c1+Ktens)*tmpcalc + rep_prssw = (c1-Ktens)*tmpcalc*Deltasw + etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot + +! else + +! endif + + end subroutine viscous_coeffs_and_rep_pressure + !======================================================================= ! Load velocity components into array for boundary updates 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/doc/source/cice_index.rst b/doc/source/cice_index.rst index 19428c87e..12bc8d32e 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -202,6 +202,7 @@ 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" @@ -515,7 +516,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 +540,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 +735,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/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index d4e209d8a..6d7d32976 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -385,7 +385,7 @@ 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`_. +Documentation `_. .. _stress-vp: @@ -405,7 +405,7 @@ An elliptical yield curve is used, with the viscosities given by \zeta = {P(1+k_t)\over 2\Delta}, .. math:: - \eta = {P(1+k_t)\over {2\Delta e^2}}, + \eta = e^{-2} \zeta, where @@ -455,29 +455,28 @@ 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(1-k_t)\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(1-k_t)\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`. @@ -721,11 +720,10 @@ Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite .. 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,\\ + + {P_R^k(1-k_t)} &=& 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}} &=& - {P(1+k_t)\over 2e^2\Delta^k}D_S^k,\end{aligned} + \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`. diff --git a/doc/source/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index ce5c2ef41..955328bf5 100644 --- a/doc/source/user_guide/ug_testing.rst +++ b/doc/source/user_guide/ug_testing.rst @@ -1121,7 +1121,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 +1130,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 From 9f72e876107ccb18c2e3a5c09fe4a50caf8e8250 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Thu, 14 Oct 2021 10:08:20 -0700 Subject: [PATCH 09/19] Deprecate gx1 CORE forcing option (#643) * deprecate gx1 core forcing option --- cicecore/cicedynB/general/ice_forcing.F90 | 15 +++++++++++---- configuration/scripts/options/set_nml.gx1coreii | 10 ---------- configuration/scripts/tests/base_suite.ts | 1 - doc/source/developer_guide/dg_forcing.rst | 15 --------------- doc/source/user_guide/ug_case_settings.rst | 1 - 5 files changed, 11 insertions(+), 31 deletions(-) delete mode 100644 configuration/scripts/options/set_nml.gx1coreii 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/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/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 2c75bffe1..3dc4905b3 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -70,7 +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/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/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 2b4e26c72..23e6951fc 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -542,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", "" From 2207d8843c3497b2d1c4f72d4bce90f0f8d11cef Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Wed, 20 Oct 2021 08:23:37 -0700 Subject: [PATCH 10/19] Update test features: add create_fails script, QC process test, and wait on baseline results (#644) * add create_fails.csh and update documentation * update documentation * update qcnonbfb settings * add qc test capability and update bfbcomp logic to add check that the baseline results has completed * add PEND for bfbcomp result * update regression testing to use QC (instead of log compare) for qcchk test cases * update create_fails.csh * update create_fails log output --- cice.setup | 7 ++- .../scripts/machines/env.cheyenne_gnu | 7 +++ .../scripts/machines/env.cheyenne_intel | 7 +++ .../scripts/machines/env.cheyenne_pgi | 7 +++ .../scripts/options/set_nml.qc_nonbfb | 8 ---- .../scripts/options/set_nml.qcnonbfb | 16 +++++++ configuration/scripts/tests/QC/cice.t-test.py | 10 +++- .../scripts/tests/QC/gen_qc_cases.csh | 4 +- configuration/scripts/tests/baseline.script | 48 +++++++++++++++++++ configuration/scripts/tests/create_fails.csh | 24 ++++++++++ configuration/scripts/tests/decomp_suite.ts | 2 - configuration/scripts/tests/first_suite.ts | 1 + configuration/scripts/tests/prod_suite.ts | 8 ++-- configuration/scripts/tests/reprosum_suite.ts | 1 - configuration/scripts/tests/test_qcchk.script | 36 ++++++++++++++ .../scripts/tests/test_qcchkf.script | 36 ++++++++++++++ doc/source/developer_guide/dg_scripts.rst | 11 +++-- doc/source/user_guide/ug_testing.rst | 11 ++++- 18 files changed, 218 insertions(+), 26 deletions(-) delete mode 100644 configuration/scripts/options/set_nml.qc_nonbfb create mode 100644 configuration/scripts/options/set_nml.qcnonbfb create mode 100755 configuration/scripts/tests/create_fails.csh create mode 100644 configuration/scripts/tests/test_qcchk.script create mode 100644 configuration/scripts/tests/test_qcchkf.script diff --git a/cice.setup b/cice.setup index 3cae17d74..be9266dd2 100755 --- a/cice.setup +++ b/cice.setup @@ -480,6 +480,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 @@ -919,7 +920,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 #------------------------------------------------------------ @@ -1144,7 +1145,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 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_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/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/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/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/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/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index 955328bf5..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 ~~~~~~~~~~~~~~~~~~~~~~~~ From 4147db49dbc93ba93cbb5560b1c299f96a98f4d7 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Thu, 21 Oct 2021 11:14:25 -0400 Subject: [PATCH 11/19] ice_dyn_vp: express rheology in terms of bulk and shear viscosities (#647) * Added array for eta and changed name of zetaD...BFB * Added array for replacement pressure...BFB * Replacement pressure now used on calc...BFB * eta is now used in the calc...BFB * Ktens included in zeta,eta and rep_prs...roundoff errors...BFB only if Ktens=0 * Small modif to viscous coeff subroutine for passing logical capping * Modifs to viscous_coeffs suroutine to be used for VP and EVP * Further modifs to viscous_coeff subroutine * cosmetic change: order of calc is ne-nw-sw-se * vp solver also uses viscous_coeffs_and_rep_pressure subroutine * Minor modifications following PR review --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 12 +- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 57 +++-- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 242 ++++++++++-------- 3 files changed, 175 insertions(+), 136 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 775f2caf1..8f3fc4910 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -656,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)' !----------------------------------------------------------------- @@ -663,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) @@ -694,13 +697,14 @@ subroutine stress (nx_block, ny_block, & call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),& Deltane, Deltanw, & - Deltase, Deltasw, & + Deltasw, Deltase, & zetax2ne, zetax2nw, & - zetax2se, zetax2sw, & + zetax2sw, zetax2se, & etax2ne, etax2nw, & - etax2se, etax2sw, & + etax2sw, etax2se, & rep_prsne, rep_prsnw, & - rep_prsse, rep_prssw ) + rep_prssw, rep_prsse, & + capping) !----------------------------------------------------------------- ! the stresses ! kg/s^2 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 167bc6a2c..5e14d9686 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -1375,53 +1375,64 @@ end subroutine strain_rates subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & Deltane, Deltanw, & - Deltase, Deltasw, & + Deltasw, Deltase, & zetax2ne, zetax2nw, & - zetax2se, zetax2sw, & + zetax2sw, zetax2se, & etax2ne, etax2nw, & - etax2se, etax2sw, & + etax2sw, etax2se, & rep_prsne, rep_prsnw,& - rep_prsse, rep_prssw ) + 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, Deltase, Deltasw ! Delta at each corner + Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner + logical , intent(in):: capping + real (kind=dbl_kind), intent(out):: & - zetax2ne, zetax2nw, zetax2se, zetax2sw, & ! zetax2 at each corner - etax2ne, etax2nw, etax2se, etax2sw, & ! etax2 at each corner - rep_prsne, rep_prsnw, rep_prsse, rep_prssw ! replacement pressure + 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) :: & - tmpcalc + 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 - tmpcalc = strength/max(Deltane,tinyarea) ! northeast - zetax2ne = (c1+Ktens)*tmpcalc - rep_prsne = (c1-Ktens)*tmpcalc*Deltane + 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 = ecci*zetax2ne ! CHANGE FOR e_plasticpot - tmpcalc = strength/max(Deltanw,tinyarea) ! northwest - zetax2nw = (c1+Ktens)*tmpcalc - rep_prsnw = (c1-Ktens)*tmpcalc*Deltanw + zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest + rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot - tmpcalc = strength/max(Deltase,tinyarea) ! southeast - zetax2se = (c1+Ktens)*tmpcalc - rep_prsse = (c1-Ktens)*tmpcalc*Deltase - etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot - - tmpcalc = strength/max(Deltasw,tinyarea) ! southwest - zetax2sw = (c1+Ktens)*tmpcalc - rep_prssw = (c1-Ktens)*tmpcalc*Deltasw + zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest + rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot + zetax2se = (c1+Ktens)*tmpcalcse ! southeast + rep_prsse = (c1-Ktens)*tmpcalcse*Deltase + etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot + ! else ! endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 860865dba..1a6c68548 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -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 @@ -1312,7 +1327,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 +1357,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 +1406,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) + + 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) = 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 + 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 +1552,7 @@ subroutine matvec (nx_block, ny_block, & cxm , cym , & uvel , vvel , & vrel , Cb , & - zetaD , & + zetax2 , etax2 , & umassdti, fm , & uarear , & Au , Av) @@ -1572,7 +1590,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 +1666,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 +2010,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 +2032,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 +2221,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 +2678,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 +2694,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 +2806,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 +2877,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 +2913,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 +3079,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 +3090,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 +3199,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 +3271,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 +3295,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 +3408,7 @@ subroutine pgmres (zetaD , & end do ! Call preconditioner - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & workspace_x , workspace_y, & @@ -3451,7 +3474,7 @@ end subroutine pgmres ! ! authors: Philippe Blain, ECCC - subroutine precondition(zetaD , & + subroutine precondition(zetax2 , etax2, & Cb , vrel , & umassdti , & vx , vy , & @@ -3460,8 +3483,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 +3549,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 , & From b626c49bdd45f571934916501e46f0f10678324f Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Wed, 3 Nov 2021 16:47:12 -0400 Subject: [PATCH 12/19] doc: user_guide: add missing word in Troubleshooting / Restarts (#652) --- doc/source/user_guide/ug_troubleshooting.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/user_guide/ug_troubleshooting.rst b/doc/source/user_guide/ug_troubleshooting.rst index cd8f1acaf..dc97be63f 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 From 3ea1a50313d560cff3fd41e896cec9a97aa8c238 Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Fri, 5 Nov 2021 15:53:36 -0400 Subject: [PATCH 13/19] doc: remove mention of renamed variable 'forcing_diag' (#657) 'forcing_diag' was renamed to 'debug_forcing' in d6eb125 (Add new unit tests sumchk and bcstchk and update tests (#606), 2021-06-09), but this instance was not renamed in the doc. --- doc/source/user_guide/ug_troubleshooting.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/user_guide/ug_troubleshooting.rst b/doc/source/user_guide/ug_troubleshooting.rst index dc97be63f..d20b14ffc 100644 --- a/doc/source/user_guide/ug_troubleshooting.rst +++ b/doc/source/user_guide/ug_troubleshooting.rst @@ -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**) From aaf392b9094af918a38136aed62186b5848c1c1d Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Tue, 9 Nov 2021 10:44:34 -0500 Subject: [PATCH 14/19] Implementation of plastic potential for VP and EVP. (#649) * Implementation of plastic potential * Modified ice_in * Corrected minor issue for variable declaration * Modifs to the doc for plastic potential * Minor modifs to the doc * More details in doc about the plastic potential * Added journal TC=the cryosphere in master_list.bib * Modified the doc following comments from Elizabeth and Philippe --- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 32 +-- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 5 +- cicecore/cicedynB/general/ice_init.F90 | 32 +-- configuration/scripts/ice_in | 3 +- configuration/scripts/options/set_nml.alt03 | 2 +- doc/source/cice_index.rst | 4 +- doc/source/master_list.bib | 13 +- doc/source/science_guide/sg_dynamics.rst | 198 +++++++++--------- 8 files changed, 159 insertions(+), 130 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 5e14d9686..23251b2d1 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -63,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 @@ -220,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)' @@ -230,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 @@ -1352,10 +1354,10 @@ 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 @@ -1419,19 +1421,19 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & zetax2ne = (c1+Ktens)*tmpcalcne ! northeast rep_prsne = (c1-Ktens)*tmpcalcne*Deltane - etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot + etax2ne = epp2i*zetax2ne zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw - etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot + etax2nw = epp2i*zetax2nw zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw - etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot + etax2sw = epp2i*zetax2sw zetax2se = (c1+Ktens)*tmpcalcse ! southeast rep_prsse = (c1-Ktens)*tmpcalcse*Deltase - etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot + etax2se = epp2i*zetax2se ! else diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 1a6c68548..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 @@ -1320,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 , & diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 2e67af51c..b299ef77f 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 @@ -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) @@ -1440,9 +1443,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 diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 443ff1cbb..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 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/doc/source/cice_index.rst b/doc/source/cice_index.rst index 12bc8d32e..85acbece3 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -205,7 +205,9 @@ either Celsius or Kelvin units). "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" 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 6d7d32976..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 = e^{-2} \zeta, + \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``). + +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`. -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. +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. -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``. +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,7 +463,7 @@ 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} &=& {\zeta \over T} D_D, \\ + + {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} &=& @@ -466,7 +474,7 @@ 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} &=& {\zeta^k\over T} D_D^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} &=& @@ -480,6 +488,78 @@ the viscosity terms in the subcycling. Choices of the parameters used to define :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 @@ -667,77 +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)} &=& 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`. - -.. _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. From 47c2d97f3bcdb203c65206b8bad90c398fa61fb4 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Mon, 15 Nov 2021 18:07:09 -0800 Subject: [PATCH 15/19] update icepack, rename snwITDrdg to snwitdrdg (#658) --- .../scripts/options/{set_nml.snwITDrdg => set_nml.snwitdrdg} | 0 configuration/scripts/tests/base_suite.ts | 4 ++-- icepack | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) rename configuration/scripts/options/{set_nml.snwITDrdg => set_nml.snwitdrdg} (100%) 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/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 3dc4905b3..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 diff --git a/icepack b/icepack index f9c9e480f..152bd701e 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit f9c9e480f6ce482317734be80719178c8e1b5121 +Subproject commit 152bd701e0cf3ec4385e5ce81918ba94e7a791cb From 2ccb8f2fcf51bcb3b764ff26becb9990676a4b9d Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Fri, 19 Nov 2021 18:51:53 -0700 Subject: [PATCH 16/19] Change max_blocks for rake tests on izumi (nothread). (#665) * Fix some raketests for izumi * fix some rake tests --- configuration/scripts/tests/nothread_suite.ts | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 410f39fffed0fc44e9ff319ecca46a157ad6a5a4 Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Fri, 26 Nov 2021 12:01:54 -0500 Subject: [PATCH 17/19] Makefile: make Fortran object files depend on their dependency files (#667) When 'make' is invoked on the CICE Makefile, the first thing it does is to try to make the included dependency files (*.d) (which are in fact Makefiles themselves) [1], in alphabetical order. The rule to make the dep files have the dependency generator, 'makdep', as a prerequisite, so when processing the first dep file, make notices 'makdep' does not exist and proceeds to build it. If for whatever reason this compilation fails, make will then proceed to the second dep file, notice that it recently tried and failed to build its dependency 'makdep', give up on the second dep file, proceed to the third, and so on. In the end, no dep file is produced. Make then restarts itself and proceeds to build the code, which of course fails catastrophically because the Fortran source files are not compiled in the right order because the dependency files are missing. To avoid that, add a dependency on the dep file to the rules that make the object file out of the Fortran source files. Since old-fashioned suffix rules cannot have their own prerequisites [2], migrate the rules for the Fortran source files to use pattern rules [3] instead. While at it, also migrate the rule for the C source files. With this new dependency, the builds abort early, before trying to compile the Fortran sources, making it easier to understand what has gone wrong. Since we do not use suffix rules anymore, remove the '.SUFFIXES' line that indicates which extension to use suffix rules for (but keep the line that eliminates all default suffix rules). [1] https://www.gnu.org/software/make/manual/html_node/Remaking-Makefiles.html [2] https://www.gnu.org/software/make/manual/html_node/Suffix-Rules.html [3] https://www.gnu.org/software/make/manual/html_node/Pattern-Rules.html#Pattern-Rules --- configuration/scripts/Makefile | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) 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: From 162aee96bc344ed5759901b010876960a5da57aa Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Fri, 26 Nov 2021 10:12:39 -0800 Subject: [PATCH 18/19] Fix multi-pe advection=none bug (#664) * update parsing scripts to improve robustness, fix multi-pe advection=none * Update cice script to improve performance including minor refactoring of parse_namelist and parse_settings to reduce cost and ability to use already setup ice_in file from a prior case in the suite. Added commented out timing ability in cice.setup. Change test default to PEND from FAIL. * fix cice.setup for case * add sedbak implementation to support Mac sed * s/spend/spent --- cice.setup | 73 ++++++++++++++++--- cicecore/cicedynB/general/ice_init.F90 | 7 +- configuration/scripts/options/set_env.box2001 | 1 - configuration/scripts/parse_namelist.sh | 18 ++--- .../scripts/parse_namelist_from_env.sh | 3 +- configuration/scripts/parse_settings.sh | 34 +++++---- 6 files changed, 98 insertions(+), 38 deletions(-) delete mode 100644 configuration/scripts/options/set_env.box2001 diff --git a/cice.setup b/cice.setup index be9266dd2..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 @@ -824,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 @@ -934,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} @@ -965,7 +979,6 @@ EOF1 cat >! ${fsmods} << EOF1 # cice.setup settings - setenv ICE_SANDBOX ${ICE_SANDBOX} setenv ICE_SCRIPTS ${ICE_SCRIPTS} setenv ICE_CASENAME ${casename} @@ -1034,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" @@ -1078,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 @@ -1162,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 @@ -1176,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 @@ -1218,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/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index b299ef77f..7485cbe23 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -958,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" @@ -1467,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 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/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 From 2b851261b59a9fae63a96b4e52aea3b96a0ec5d9 Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Mon, 29 Nov 2021 18:12:10 -0500 Subject: [PATCH 19/19] nuopc/cmeps driver updates (#668) * add debug_model feature * add required variables and calls for tr_snow --- cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 | 49 +++++++++--- cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 | 80 +++++++++++++++++-- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 15 ++++ .../drivers/nuopc/cmeps/ice_import_export.F90 | 2 +- 4 files changed, 129 insertions(+), 17 deletions(-) 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/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index a832e7bdf..cad480dd9 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -88,6 +88,7 @@ module ice_comp_nuopc integer :: nthrds ! Number of threads to use in this component integer :: dbug = 0 + logical :: profile_memory = .false. integer , parameter :: debug_import = 0 ! internal debug level integer , parameter :: debug_export = 0 ! internal debug level character(*), parameter :: modName = "(ice_comp_nuopc)" @@ -157,6 +158,10 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) type(ESMF_State) :: importState, exportState type(ESMF_Clock) :: clock integer, intent(out) :: rc + + logical :: isPresent, isSet + character(len=64) :: value + character(len=char_len_long) :: logmsg !-------------------------------- rc = ESMF_SUCCESS @@ -166,6 +171,14 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) acceptStringList=(/"IPDv01p"/), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + profile_memory = .false. + call NUOPC_CompAttributeGet(gcomp, name="ProfileMemory", value=value, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) profile_memory=(trim(value)=="true") + write(logmsg,*) profile_memory + call ESMF_LogWrite('CICE_cap:ProfileMemory = '//trim(logmsg), ESMF_LOGMSG_INFO) + end subroutine InitializeP0 !=============================================================================== @@ -1049,7 +1062,9 @@ subroutine ModelAdvance(gcomp, rc) ! Advance cice and timestep update !-------------------------------- + if(profile_memory) call ESMF_VMLogMemInfo("Entering CICE_Run : ") call CICE_Run() + if(profile_memory) call ESMF_VMLogMemInfo("Leaving CICE_Run : ") !-------------------------------- ! Create export state diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 62ff2727d..f8627d690 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -582,7 +582,7 @@ subroutine ice_import( importState, rc ) rhoa(i,j,iblk) = inst_pres_height_lowest / & (287.058_ESMF_KIND_R8*(1._ESMF_KIND_R8+0.608_ESMF_KIND_R8*Qa(i,j,iblk))*Tair(i,j,iblk)) else - rhoa(i,j,iblk) = 0._ESMF_KIND_R8 + rhoa(i,j,iblk) = 1.2_ESMF_KIND_R8 endif end do !i end do !j