diff --git a/.github/workflows/test-cice.yml b/.github/workflows/test-cice.yml new file mode 100644 index 000000000..1fdd8188d --- /dev/null +++ b/.github/workflows/test-cice.yml @@ -0,0 +1,135 @@ +name: GHActions + +# This workflow is triggered on pushes, pull-requeust, and releases +# ghactions* branch names will trigger this to support development testing +# To Do: get it working with bash and ubuntu + +on: + push: + branches: + - master + - 'CICE*' + - 'ghactions*' + pull_request: + release: + types: + - created + +defaults: + run: + shell: /bin/csh {0} + +jobs: + build: + name: "CICETesting" + runs-on: ${{ matrix.os }} + strategy: + matrix: +# os: [macos-latest, ubuntu-latest] + os: [macos-latest] +# os: [ubuntu-latest] + include: + - os: macos-latest + envdef: macos + minicond: Miniconda3-latest-MacOSX-x86_64.sh +# - os: ubuntu-latest +# envdef: linux +# minicond: Miniconda3-latest-Linux-x86_64.sh + steps: + - name: reset macos toolchain to commandlinetools + shell: /bin/bash {0} + if: contains( matrix.envdef, 'macos') + run: | + sudo xcode-select -r + sudo xcode-select -s /Library/Developer/CommandLineTools + echo "xcrun --show-sdk-path: $(xcrun --show-sdk-path)" + echo "xcode-select -p: $(xcode-select -p)" + - name: system info + shell: /bin/bash {0} + run: | + type wget + type curl + type csh + echo "readlink \$(which csh): $(python -c 'import os, sys; print os.path.realpath(sys.argv[1])' $(which csh))" + echo "csh --version: $(csh --version)" + echo "uname -a: $(uname -a)" + echo "sw_vers: $(sw_vers)" + echo "HOME: $HOME" + echo "GITHUB_WORKSPACE: $GITHUB_WORKSPACE" + echo "OS: ${{ matrix.os }}" + echo "ENVDEF: ${{ matrix.envdef }}" + echo "MINICOND: ${{ matrix.minicond }}" + - name : install miniconda + shell: /bin/bash {0} + run: | + wget https://repo.anaconda.com/miniconda/${{ matrix.minicond }} -O ~/miniconda.sh + bash ~/miniconda.sh -b -p $HOME/miniconda + - name: clone + uses: actions/checkout@v2 + with: + submodules: 'recursive' + - name: link + run: | + ln -s ${GITHUB_WORKSPACE}/../CICE ${HOME}/cice +# ls -al ${HOME}/ +# ls -al ${GITHUB_WORKSPACE}/ + - name: setup conda env + shell: /bin/bash {0} + run: | + cd $HOME && mkdir -p cice-dirs/runs cice-dirs/baseline cice-dirs/input + source $HOME/miniconda/bin/activate + conda init tcsh + cd $HOME/cice + conda env create -f configuration/scripts/machines/environment.yml + - name: check conda env + run: | + conda activate cice && which mpicc && which mpifort && which make + mpifort --version + mpicc --version + make --version + - name: check setup case + run: | + cd $HOME/cice + ./cice.setup -m conda -e ${{ matrix.envdef }} -c case0 --pes 1x1 -s diag1 + - name: check setup test + run: | + cd $HOME/cice + ./cice.setup -m conda -e ${{ matrix.envdef }} --test smoke --testid c0 +# - name: compile case +# run: | +# cd $HOME/cice +# ./cice.setup -m conda -e ${{ matrix.envdef }} -c case1 +# cd case1 +# ./cice.build + - name: download input data + run: | + cd $HOME/cice-dirs/input + wget https://zenodo.org/record/3728358/files/CICE_data_gx3_grid_ic-20200320.tar.gz && tar xvfz CICE_data_gx3_grid_ic-20200320.tar.gz + wget https://zenodo.org/record/3728362/files/CICE_data_gx3_forcing_NCAR_bulk-20200320.tar.gz && tar xvfz CICE_data_gx3_forcing_NCAR_bulk-20200320.tar.gz + wget https://zenodo.org/record/3728364/files/CICE_data_gx3_forcing_JRA55-20200320.tar.gz && tar xvfz CICE_data_gx3_forcing_JRA55-20200320.tar.gz + pwd + ls -alR +# - name: run case +# run: | +# cd $HOME/cice +# cd case1 +# ./cice.run + - name: run suite + run: | + cd $HOME/cice + ./cice.setup -m conda -e ${{ matrix.envdef }} --suite travis_suite --testid ${{ matrix.os }} + - name: write output + run: | + cd $HOME/cice + ./.github/workflows/write_logfiles.csh + cd testsuite.${{ matrix.os }} + ./results.csh + - name: successful run + if: ${{ success() }} + run: | + echo "${{ job.name }} PASSED" + - name: trap failure + if: ${{ failure() }} + run: | + echo "${{ job.name }} FAILED" + exit 99 diff --git a/.github/workflows/write_logfiles.csh b/.github/workflows/write_logfiles.csh new file mode 100755 index 000000000..a6777dec6 --- /dev/null +++ b/.github/workflows/write_logfiles.csh @@ -0,0 +1,10 @@ +#!/bin/csh + +#echo "hello" + +foreach logfile (case*/logs/cice.runlog* testsuite.*/*/logs/cice.runlog*) + echo "### ${logfile} ###" + tail -20 $logfile + echo " " +end + diff --git a/.gitmodules b/.gitmodules index 8a773d230..b84a13b43 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,4 @@ [submodule "icepack"] path = icepack - url = https://github.com/NOAA-EMC/Icepack + #url = https://github.com/NOAA-EMC/Icepack + url = https://github.com/DeniseWorthen/Icepack diff --git a/.travis.yml b/.travis.yml index f8f5aeadc..994079e98 100644 --- a/.travis.yml +++ b/.travis.yml @@ -25,7 +25,11 @@ install: # Fetch forcing data - "wget https://zenodo.org/record/3728362/files/CICE_data_gx3_forcing_NCAR_bulk-20200320.tar.gz && - tar xvfz CICE_data_gx3_forcing_NCAR_bulk-20200320.tar.gz -C ~" + tar xvfz CICE_data_gx3_forcing_NCAR_bulk-20200320.tar.gz -C ~" + + # Fetch jra55_gx3 forcing data + - "wget https://zenodo.org/record/3728364/files/CICE_data_gx3_forcing_JRA55-20200320.tar.gz && + tar xvfz CICE_data_gx3_forcing_JRA55-20200320.tar.gz -C ~" # Mirror entire data folder #- "lftp ftp://anonymous:travis@travis-ci.org@ftp.cgd.ucar.edu diff --git a/LICENSE.pdf b/LICENSE.pdf index da37344cf..5d6b29280 100644 Binary files a/LICENSE.pdf and b/LICENSE.pdf differ diff --git a/README.md b/README.md index a584e8ac9..457714b8a 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ -[![Build Status](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=master)](https://travis-ci.org/CICE-Consortium/CICE) +[![GHActions](https://github.com/CICE-Consortium/CICE/workflows/GHActions/badge.svg)](https://github.com/CICE-Consortium/CICE/actions) [![Documentation Status](https://readthedocs.org/projects/cice-consortium-cice/badge/?version=master)](http://cice-consortium-cice.readthedocs.io/en/master/?badge=master) [![lcov](https://img.shields.io/endpoint?url=https://apcraig.github.io/coverage.json)](https://apcraig.github.io) diff --git a/cicecore/cicedynB/analysis/ice_diagnostics.F90 b/cicecore/cicedynB/analysis/ice_diagnostics.F90 index 40da6cb64..cff544cd4 100644 --- a/cicecore/cicedynB/analysis/ice_diagnostics.F90 +++ b/cicecore/cicedynB/analysis/ice_diagnostics.F90 @@ -1561,11 +1561,11 @@ subroutine print_state(plabel,i,j,iblk) real (kind=dbl_kind) :: & eidebug, esdebug, & - qi, qs, Tsnow, & + qi, qs, Tsnow, si, & rad_to_deg, puny, rhoi, lfresh, rhos, cp_ice integer (kind=int_kind) :: n, k, nt_Tsfc, nt_qice, nt_qsno, nt_fsd, & - nt_isosno, nt_isoice + nt_isosno, nt_isoice, nt_sice logical (kind=log_kind) :: tr_fsd, tr_iso @@ -1576,7 +1576,7 @@ subroutine print_state(plabel,i,j,iblk) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, & - nt_qsno_out=nt_qsno, nt_fsd_out=nt_fsd, & + nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, & nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) call icepack_query_parameters( & rad_to_deg_out=rad_to_deg, puny_out=puny, rhoi_out=rhoi, lfresh_out=lfresh, & @@ -1621,7 +1621,6 @@ subroutine print_state(plabel,i,j,iblk) enddo ! n - eidebug = c0 do n = 1,ncat do k = 1,nilyr @@ -1654,6 +1653,14 @@ subroutine print_state(plabel,i,j,iblk) write(nu_diag,*) 'qsnow(i,j)',esdebug write(nu_diag,*) ' ' + do n = 1,ncat + do k = 1,nilyr + si = trcrn(i,j,nt_sice+k-1,n,iblk) + write(nu_diag,*) 'sice, cat ',n,' layer ',k, si + enddo + enddo + write(nu_diag,*) ' ' + write(nu_diag,*) 'uvel(i,j)',uvel(i,j,iblk) write(nu_diag,*) 'vvel(i,j)',vvel(i,j,iblk) diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index 660676a64..1aa2515a4 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/cicecore/cicedynB/analysis/ice_history.F90 @@ -925,12 +925,12 @@ subroutine init_hist (dt) ns1, f_strinty) call define_hist_field(n_taubx,"taubx","N/m^2",ustr2D, ucstr, & - "basal (seabed) stress (x)", & + "seabed (basal) stress (x)", & "positive is x direction on U grid", c1, c0, & ns1, f_taubx) call define_hist_field(n_tauby,"tauby","N/m^2",ustr2D, ucstr, & - "basal (seabed) stress (y)", & + "seabed (basal) stress (y)", & "positive is y direction on U grid", c1, c0, & ns1, f_tauby) diff --git a/cicecore/cicedynB/analysis/ice_history_bgc.F90 b/cicecore/cicedynB/analysis/ice_history_bgc.F90 index 1ae572b30..67b23904e 100644 --- a/cicecore/cicedynB/analysis/ice_history_bgc.F90 +++ b/cicecore/cicedynB/analysis/ice_history_bgc.F90 @@ -267,7 +267,7 @@ module ice_history_bgc subroutine init_hist_bgc_2D use ice_broadcast, only: broadcast_scalar - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_communicate, only: my_task, master_task use ice_history_shared, only: tstr2D, tcstr, define_hist_field, & f_fsalt, f_fsalt_ai, f_sice @@ -781,6 +781,7 @@ subroutine init_hist_bgc_2D if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_iso(1:1) /= 'x') then do n=1,n_iso @@ -1780,6 +1781,7 @@ subroutine init_hist_bgc_2D "distance from ice bottom to brine surface", c1, c0, & ns, f_hbri) + endif ! histfreq(ns) /= 'x' enddo ! nstreams endif ! tr_aero, etc @@ -1790,7 +1792,7 @@ end subroutine init_hist_bgc_2D subroutine init_hist_bgc_3Dc - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_history_shared, only: tstr3Dc, tcstr, define_hist_field integer (kind=int_kind) :: ns @@ -1802,18 +1804,19 @@ subroutine init_hist_bgc_3Dc if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (tr_brine) then + if (tr_brine) then ! 3D (category) variables must be looped separately do ns = 1, nstreams - if (f_fbri(1:1) /= 'x') & - call define_hist_field(n_fbri,"fbrine","1",tstr3Dc, tcstr, & + if (histfreq(ns) /= 'x') then + if (f_fbri(1:1) /= 'x') & + call define_hist_field(n_fbri,"fbrine","1",tstr3Dc, tcstr, & "brine tracer fraction of ice volume, cat", & - "none", c1, c0, & - ns, f_fbri) + "none", c1, c0, ns, f_fbri) + endif ! histfreq /= 'x' enddo ! ns - endif + endif ! tr_brine end subroutine init_hist_bgc_3Dc @@ -1821,7 +1824,7 @@ end subroutine init_hist_bgc_3Dc subroutine init_hist_bgc_3Db - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams,histfreq use ice_history_shared, only: tstr3Db, tcstr, define_hist_field integer (kind=int_kind) :: ns @@ -1841,6 +1844,7 @@ subroutine init_hist_bgc_3Db if (z_tracers .or. solve_zsal) then do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_bTin(1:1) /= 'x') & call define_hist_field(n_bTin,"bTizn","C",tstr3Db, tcstr, & @@ -1873,6 +1877,7 @@ subroutine init_hist_bgc_3Db "internal ice PAR", "on bio interface grid", c1, c0, & ns, f_zfswin) + endif ! histfreq(ns) /= 'x' enddo ! ns endif ! z_tracers or solve_zsal @@ -2012,9 +2017,10 @@ subroutine accum_hist_bgc (iblk) ! increment field !--------------------------------------------------------------- - if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then - ! 2d bgc fields + ! 2d bgc fields + if (allocated(a2D)) then + if (tr_iso .or. tr_aero .or. tr_brine .or. solve_zsal .or. skl_bgc) then ! zsalinity if (f_fzsal (1:1) /= 'x') & @@ -2635,11 +2641,12 @@ subroutine accum_hist_bgc (iblk) call accum_hist_field(n_hbri, iblk, & hbri(:,:,iblk), a2D) - endif ! 2d bgc tracers, tr_aero, tr_brine, solve_zsal, skl_bgc - + 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 ! 3Dc bgc category fields @@ -2647,7 +2654,9 @@ subroutine accum_hist_bgc (iblk) call accum_hist_field(n_fbri-n2D, iblk, ncat_hist, & trcrn(:,:,nt_fbri,1:ncat_hist,iblk), a3Dc) endif + endif ! allocated(a3Dc) + if (allocated(a3Db)) then if (z_tracers .or. solve_zsal) then ! 3Db category fields @@ -2754,8 +2763,10 @@ subroutine accum_hist_bgc (iblk) workz(:,:,1:nzblyr), a3Db) endif - endif ! 3Db fields + endif ! 3Db fields + endif ! allocated(a3Db) + if (allocated(a3Da)) then if (z_tracers) then ! 3Da category fields @@ -3189,7 +3200,8 @@ subroutine accum_hist_bgc (iblk) workz2(:,:,1:nzalyr), a3Da) endif - endif ! z_tracers, 3Da tracers + endif ! z_tracers, 3Da tracers + endif ! allocated(a3Da) end subroutine accum_hist_bgc @@ -3197,7 +3209,7 @@ end subroutine accum_hist_bgc subroutine init_hist_bgc_3Da - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_history_shared, only: tstr3Da, tcstr, define_hist_field integer (kind=int_kind) :: ns, n @@ -3216,6 +3228,7 @@ subroutine init_hist_bgc_3Da if (z_tracers) then do ns = 1, nstreams + if (histfreq(ns) /= 'x') then !---------------------------------------------------------------------------- ! snow+bio grid ==> @@ -3439,6 +3452,7 @@ subroutine init_hist_bgc_3Da "other bulk nitrogen pool in cat 1", "snow+bio grid", c1, c0, & ns, f_bgc_PON_cat1) + endif ! histfreq(ns) /= 'x' enddo !ns endif ! z_tracers diff --git a/cicecore/cicedynB/analysis/ice_history_mechred.F90 b/cicecore/cicedynB/analysis/ice_history_mechred.F90 index 78edd7a31..a20df5fb0 100644 --- a/cicecore/cicedynB/analysis/ice_history_mechred.F90 +++ b/cicecore/cicedynB/analysis/ice_history_mechred.F90 @@ -81,7 +81,7 @@ module ice_history_mechred subroutine init_hist_mechred_2D use ice_broadcast, only: broadcast_scalar - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_communicate, only: my_task, master_task use ice_history_shared, only: tstr2D, tcstr, define_hist_field @@ -157,6 +157,7 @@ subroutine init_hist_mechred_2D ! 2D variables do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_alvl(1:1) /= 'x') & call define_hist_field(n_alvl,"alvl","1",tstr2D, tcstr, & @@ -203,6 +204,7 @@ subroutine init_hist_mechred_2D "none", secday*c100, c0, & ns, f_opening) + endif ! histfreq(ns) /= 'x' enddo ! nstreams end subroutine init_hist_mechred_2D @@ -211,7 +213,7 @@ end subroutine init_hist_mechred_2D subroutine init_hist_mechred_3Dc - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_history_shared, only: tstr3Dc, tcstr, define_hist_field integer (kind=int_kind) :: ns @@ -228,6 +230,7 @@ subroutine init_hist_mechred_3Dc file=__FILE__, line=__LINE__) do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_ardgn(1:1) /= 'x') & call define_hist_field(n_ardgn,"ardgn","1",tstr3Dc, tcstr, & @@ -295,6 +298,7 @@ subroutine init_hist_mechred_3Dc "none", c1, c0, & ns, f_vraftn) + endif ! histfreq(ns) /= 'x' enddo ! ns end subroutine init_hist_mechred_3Dc @@ -332,6 +336,7 @@ subroutine accum_hist_mechred (iblk) !--------------------------------------------------------------- ! 2D fields + if (allocated(a2D)) then if (f_alvl(1:1)/= 'x') & call accum_hist_field(n_alvl, iblk, & @@ -354,7 +359,10 @@ subroutine accum_hist_mechred (iblk) if (f_opening(1:1) /= 'x') & call accum_hist_field(n_opening, iblk, opening(:,:,iblk), a2D) + endif ! allocated(a2D) + ! 3D category fields + if (allocated(a3Dc)) then if (f_ardgn(1:1)/= 'x') & call accum_hist_field(n_ardgn-n2D, iblk, ncat_hist, & @@ -391,6 +399,7 @@ subroutine accum_hist_mechred (iblk) if (f_vraftn(1:1)/= 'x') & call accum_hist_field(n_vraftn-n2D, iblk, ncat_hist, & vraftn(:,:,1:ncat_hist,iblk), a3Dc) + endif ! allocated(a3Dc) end subroutine accum_hist_mechred diff --git a/cicecore/cicedynB/analysis/ice_history_pond.F90 b/cicecore/cicedynB/analysis/ice_history_pond.F90 index ebef84483..de10eb9fb 100644 --- a/cicecore/cicedynB/analysis/ice_history_pond.F90 +++ b/cicecore/cicedynB/analysis/ice_history_pond.F90 @@ -66,7 +66,7 @@ module ice_history_pond subroutine init_hist_pond_2D use ice_broadcast, only: broadcast_scalar - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_communicate, only: my_task, master_task use ice_history_shared, only: tstr2D, tcstr, define_hist_field @@ -135,6 +135,7 @@ subroutine init_hist_pond_2D ! 2D variables do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_apond(1:1) /= 'x') & call define_hist_field(n_apond,"apond","1",tstr2D, tcstr, & @@ -184,6 +185,7 @@ subroutine init_hist_pond_2D "weighted by ice area", c1, c0, & ns, f_apeff_ai) + endif ! histfreq(ns) /= 'x' enddo ! nstreams endif ! tr_pond @@ -194,7 +196,7 @@ end subroutine init_hist_pond_2D subroutine init_hist_pond_3Dc - use ice_calendar, only: nstreams + use ice_calendar, only: nstreams, histfreq use ice_history_shared, only: tstr3Dc, tcstr, define_hist_field integer (kind=int_kind) :: ns @@ -210,6 +212,7 @@ subroutine init_hist_pond_3Dc ! 3D (category) variables must be looped separately do ns = 1, nstreams + if (histfreq(ns) /= 'x') then if (f_apondn(1:1) /= 'x') & call define_hist_field(n_apondn,"apondn","1",tstr3Dc, tcstr, & @@ -227,6 +230,7 @@ subroutine init_hist_pond_3Dc "none", c1, c0, & ns, f_apeffn) + endif ! histfreq(ns) /= 'x' enddo ! ns endif ! tr_pond @@ -286,6 +290,7 @@ subroutine accum_hist_pond (iblk) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) + if (allocated(a2D)) then if (tr_pond_cesm) then if (f_apond(1:1)/= 'x') & call accum_hist_field(n_apond, iblk, & @@ -374,7 +379,10 @@ subroutine accum_hist_pond (iblk) if (f_apeff_ai(1:1) /= 'x') & call accum_hist_field(n_apeff_ai, iblk, apeff_ai(:,:,iblk), a2D) + endif ! allocated(a2D) + ! 3D category fields + if (allocated(a3Dc)) then if (f_apondn (1:1) /= 'x') & call accum_hist_field(n_apondn-n2D, iblk, ncat_hist, & trcrn(:,:,nt_apnd,1:ncat_hist,iblk), a3Dc) @@ -385,6 +393,7 @@ subroutine accum_hist_pond (iblk) call accum_hist_field(n_hpondn-n2D, iblk, ncat_hist, & trcrn(:,:,nt_apnd,1:ncat_hist,iblk) & * trcrn(:,:,nt_hpnd,1:ncat_hist,iblk), a3Dc) + endif ! allocated(a3Dc) end subroutine accum_hist_pond diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index e6bb86bff..2face07c2 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -122,7 +122,8 @@ subroutine eap (dt) use ice_dyn_shared, only: fcor_blk, ndte, dtei, & denom1, uvel_init, vvel_init, arlx1i, & dyn_prep1, dyn_prep2, stepu, dyn_finish, & - basal_stress_coeff, basalstress, & + seabed_stress_factor_LKD, seabed_stress_factor_prob, & + seabed_stress_method, seabed_stress, & stack_velocity_field, unstack_velocity_field use ice_flux, only: rdg_conv, strairxT, strairyT, & strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & @@ -383,17 +384,31 @@ subroutine eap (dt) endif !----------------------------------------------------------------- - ! basal stress coefficients (landfast ice) + ! seabed stress factor Tbu (Tbu is part of Cb coefficient) !----------------------------------------------------------------- - if (basalstress) then + if (seabed_stress) then + !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks - call basal_stress_coeff (nx_block, ny_block, & - icellu (iblk), & - indxui(:,iblk), indxuj(:,iblk), & - vice(:,:,iblk), aice(:,:,iblk), & - hwater(:,:,iblk), Tbu(:,:,iblk)) + + if ( seabed_stress_method == 'LKD' ) then + + call seabed_stress_factor_LKD (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + + elseif ( seabed_stress_method == 'probabilistic' ) then + + call seabed_stress_factor_prob (nx_block, ny_block, & + icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & + icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + aicen(:,:,:,iblk), vicen(:,:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + endif + enddo !$OMP END PARALLEL DO endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 5846cf143..d8ce42681 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -40,8 +40,9 @@ module ice_dyn_evp use ice_constants, only: c0, p027, p055, p111, p166, & p222, p25, p333, p5, c1 use ice_dyn_shared, only: stepu, dyn_prep1, dyn_prep2, dyn_finish, & - ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, & - vvel_init, basal_stress_coeff, basalstress, Ktens, revp + ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, vvel_init, & + seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & + seabed_stress, Ktens, revp use ice_fileunits, only: nu_diag use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -325,17 +326,31 @@ subroutine evp (dt) endif !----------------------------------------------------------------- - ! basal stress coefficients (landfast ice) + ! seabed stress factor Tbu (Tbu is part of Cb coefficient) !----------------------------------------------------------------- - if (basalstress) then + if (seabed_stress) then + !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - call basal_stress_coeff (nx_block, ny_block, & - icellu (iblk), & - indxui(:,iblk), indxuj(:,iblk), & - vice(:,:,iblk), aice(:,:,iblk), & - hwater(:,:,iblk), Tbu(:,:,iblk)) + do iblk = 1, nblocks + + if ( seabed_stress_method == 'LKD' ) then + + call seabed_stress_factor_LKD (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + + elseif ( seabed_stress_method == 'probabilistic' ) then + + call seabed_stress_factor_prob (nx_block, ny_block, & + icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & + icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + aicen(:,:,:,iblk), vicen(:,:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + endif + enddo !$OMP END PARALLEL DO endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 index 9fac97a89..78469cc86 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 @@ -818,7 +818,7 @@ subroutine stepu_iter(NA_len,rhow, & real (kind=dbl_kind) :: tmp_str2_nw,tmp_str3_se,tmp_str4_sw, tmp_strintx real (kind=dbl_kind) :: tmp_str6_se,tmp_str7_nw,tmp_str8_sw, tmp_strinty real (kind=dbl_kind) :: waterx,watery - real (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind ! residual velocity for basal stress (m/s) + real (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind ! residual velocity for seabed stress (m/s) character(len=*), parameter :: subname = '(stepu_iter)' !--------------------------------------- @@ -881,7 +881,7 @@ subroutine stepu_last(NA_len, rhow, & use ice_kinds_mod use ice_constants, only: c0, c1 - use ice_dyn_shared, only: brlx, revp, basalstress + use ice_dyn_shared, only: brlx, revp, seabed_stress implicit none @@ -908,7 +908,7 @@ subroutine stepu_last(NA_len, rhow, & real (kind=dbl_kind) :: tmp_str2_nw,tmp_str3_se,tmp_str4_sw real (kind=dbl_kind) :: tmp_str6_se,tmp_str7_nw,tmp_str8_sw real (kind=dbl_kind) :: waterx,watery - real (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind ! residual velocity for basal stress (m/s) + real (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind ! residual velocity for seabed stress (m/s) character(len=*), parameter :: subname = '(stepu_last)' !--------------------------------------- @@ -954,8 +954,8 @@ subroutine stepu_last(NA_len, rhow, & + umassdti(iw)*(brlx*vold + revp*vvel_init(iw)) uvel(iw) = (cca*cc1 + ccb*cc2) / ab2 vvel(iw) = (cca*cc2 - ccb*cc1) / ab2 - ! calculate basal stress component for outputs - if ( basalstress ) then + ! calculate seabed stress component for outputs + if ( seabed_stress ) then taubx(iw) = -uvel(iw)*Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) tauby(iw) = -vvel(iw)*Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index d9a0919e6..f3685ed61 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -11,8 +11,8 @@ module ice_dyn_shared use ice_kinds_mod use ice_communicate, only: my_task, master_task - use ice_constants, only: c0, c1, c2, p01, p001 - use ice_constants, only: omega, spval_dbl, p5, c4 + use ice_constants, only: c0, c1, c2, c3, c4, c6 + use ice_constants, only: omega, spval_dbl, p01, p001, p5 use ice_blocks, only: nx_block, ny_block use ice_domain_size, only: max_blocks use ice_fileunits, only: nu_diag @@ -23,7 +23,8 @@ module ice_dyn_shared implicit none private public :: init_dyn, set_evp_parameters, stepu, principal_stress, & - dyn_prep1, dyn_prep2, dyn_finish, basal_stress_coeff, & + dyn_prep1, dyn_prep2, dyn_finish, & + seabed_stress_factor_LKD, seabed_stress_factor_prob, & alloc_dyn_shared, deformations, strain_rates, & stack_velocity_field, unstack_velocity_field @@ -32,7 +33,6 @@ module ice_dyn_shared integer (kind=int_kind), public :: & kdyn , & ! type of dynamics ( -1, 0 = off, 1 = evp, 2 = eap ) kridge , & ! set to "-1" to turn off ridging - ktransport , & ! set to "-1" to turn off transport ndte ! number of subcycles: ndte=dt/dte character (len=char_len), public :: & @@ -50,8 +50,10 @@ module ice_dyn_shared ! other EVP parameters character (len=char_len), public :: & - yield_curve ! 'ellipse' ('teardrop' needs further testing) - ! + yield_curve , & ! 'ellipse' ('teardrop' needs further testing) + seabed_stress_method ! method for seabed stress calculation + ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. in prep. + real (kind=dbl_kind), parameter, public :: & eyc = 0.36_dbl_kind, & ! coefficient for calculating the parameter E @@ -82,19 +84,19 @@ module ice_dyn_shared ! ice isotropic tensile strength parameter real (kind=dbl_kind), public :: & - Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) + Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) + ! seabed (basal) stress parameters and settings logical (kind=log_kind), public :: & - basalstress ! if true, basal stress for landfast on + seabed_stress ! if true, seabed stress for landfast on - ! basal stress parameters real (kind=dbl_kind), public :: & - k1, & ! 1st free parameter for landfast parameterization - k2, & ! second free parameter (N/m^3) for landfast parametrization + k1, & ! 1st free parameter for seabed1 grounding parameterization + k2, & ! second free parameter (N/m^3) for seabed1 grounding parametrization alphab, & ! alphab=Cb factor in Lemieux et al 2015 threshold_hw, & ! max water depth for grounding ! see keel data from Amundrud et al. 2004 (JGR) - u0 = 5e-5_dbl_kind ! residual velocity for basal stress (m/s) + u0 = 5e-5_dbl_kind ! residual velocity for seabed stress (m/s) !======================================================================= @@ -447,7 +449,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & dt ! time step real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & - Tbu, & ! coefficient for basal stress (N/m^2) + Tbu, & ! seabed stress factor (N/m^2) uvel_init,& ! x-component of velocity (m/s), beginning of time step vvel_init,& ! y-component of velocity (m/s), beginning of time step umassdti, & ! mass of U-cell/dt (kg/m^2 s) @@ -469,8 +471,8 @@ subroutine dyn_prep2 (nx_block, ny_block, & strocny , & ! ice-ocean stress, y-direction strintx , & ! divergence of internal ice stress, x (N/m^2) strinty , & ! divergence of internal ice stress, y (N/m^2) - taubx , & ! basal stress, x-direction (N/m^2) - tauby ! basal stress, y-direction (N/m^2) + taubx , & ! seabed stress, x-direction (N/m^2) + tauby ! seabed stress, y-direction (N/m^2) ! local variables @@ -648,7 +650,7 @@ subroutine stepu (nx_block, ny_block, & indxuj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & - Tbu, & ! coefficient for basal stress (N/m^2) + Tbu, & ! seabed stress factor (N/m^2) uvel_init,& ! x-component of velocity (m/s), beginning of timestep vvel_init,& ! y-component of velocity (m/s), beginning of timestep aiu , & ! ice fraction on u-grid @@ -672,8 +674,8 @@ subroutine stepu (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & strintx , & ! divergence of internal ice stress, x (N/m^2) strinty , & ! divergence of internal ice stress, y (N/m^2) - taubx , & ! basal stress, x-direction (N/m^2) - tauby ! basal stress, y-direction (N/m^2) + taubx , & ! seabed stress, x-direction (N/m^2) + tauby ! seabed stress, y-direction (N/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & Cw ! ocean-ice neutral drag coefficient @@ -688,7 +690,7 @@ subroutine stepu (nx_block, ny_block, & vrel , & ! relative ice-ocean velocity cca,ccb,ab2,cc1,cc2,& ! intermediate variables taux, tauy , & ! part of ocean stress term - Cb , & ! complete basal stress coeff + Cb , & ! complete seabed (basal) stress coeff rhow ! character(len=*), parameter :: subname = '(stepu)' @@ -716,7 +718,7 @@ subroutine stepu (nx_block, ny_block, & taux = vrel*waterx(i,j) ! NOTE this is not the entire tauy = vrel*watery(i,j) ! ocn stress term - Cb = Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) ! for basal stress + Cb = Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) ! for seabed stress ! revp = 0 for classic evp, 1 for revised evp cca = (brlx + revp)*umassdti(i,j) + vrel * cosw + Cb ! kg/m^2 s @@ -739,9 +741,9 @@ subroutine stepu (nx_block, ny_block, & uvel(i,j) = (cca*cc1 + ccb*cc2) / ab2 ! m/s vvel(i,j) = (cca*cc2 - ccb*cc1) / ab2 - ! calculate basal stress component for outputs + ! calculate seabed stress component for outputs if (ksub == ndte) then ! on last subcycling iteration - if ( basalstress ) then + if ( seabed_stress ) then taubx(i,j) = -uvel(i,j)*Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) tauby(i,j) = -vvel(i,j)*Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) endif @@ -854,7 +856,10 @@ subroutine dyn_finish (nx_block, ny_block, & end subroutine dyn_finish !======================================================================= -! Computes basal stress Tbu coefficients (landfast ice) +! Computes seabed (basal) stress factor Tbu (landfast ice) based on mean +! thickness and bathymetry data. LKD refers to linear keel draft. This +! parameterization assumes that the largest keel draft varies linearly +! with the mean thickness. ! ! Lemieux, J. F., B. Tremblay, F. Dupont, M. Plante, G.C. Smith, D. Dumont (2015). ! A basal stress parameterization form modeling landfast ice, J. Geophys. Res. @@ -866,13 +871,14 @@ end subroutine dyn_finish ! ! author: JF Lemieux, Philippe Blain (ECCC) ! -! note: Tbu is a part of the Cb as defined in Lemieux et al. 2015 and 2016. -! - subroutine basal_stress_coeff (nx_block, ny_block, & - icellu, & - indxui, indxuj, & - vice, aice, & - hwater, Tbu) +! note1: Tbu is a part of the Cb as defined in Lemieux et al. 2015 and 2016. +! note2: Seabed stress (better name) was called basal stress in Lemieux et al. 2015 + + subroutine seabed_stress_factor_LKD (nx_block, ny_block, & + icellu, & + indxui, indxuj, & + vice, aice, & + hwater, Tbu) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -884,23 +890,23 @@ subroutine basal_stress_coeff (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & aice , & ! concentration of ice at tracer location - vice , & ! volume per unit area of ice at tracer location - hwater ! water depth at tracer location + vice , & ! volume per unit area of ice at tracer location (m) + hwater ! water depth at tracer location (m) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & - Tbu ! coefficient for basal stress (N/m^2) + Tbu ! seabed stress factor (N/m^2) real (kind=dbl_kind) :: & au, & ! concentration of ice at u location - hu, & ! volume per unit area of ice at u location (mean thickness) - hwu, & ! water depth at u location - hcu ! critical thickness at u location + hu, & ! volume per unit area of ice at u location (mean thickness, m) + hwu, & ! water depth at u location (m) + hcu ! critical thickness at u location (m) integer (kind=int_kind) :: & i, j, ij - character(len=*), parameter :: subname = '(basal_stress_coeff)' - + character(len=*), parameter :: subname = '(seabed1_stress_coeff)' + do ij = 1, icellu i = indxui(ij) j = indxuj(ij) @@ -917,15 +923,189 @@ subroutine basal_stress_coeff (nx_block, ny_block, & ! 1- calculate critical thickness hcu = au * hwu / k1 - ! 2- calculate basal stress factor + ! 2- calculate seabed stress factor Tbu(i,j) = k2 * max(c0,(hu - hcu)) * exp(-alphab * (c1 - au)) - + endif enddo ! ij - end subroutine basal_stress_coeff + end subroutine seabed_stress_factor_LKD + +!======================================================================= +! Computes seabed (basal) stress factor Tbu (landfast ice) based on +! probability of contact between the ITD and the seabed. The water depth +! could take into account variations of the SSH. In the simplest +! formulation, hwater is simply the value of the bathymetry. To calculate +! the probability of contact, it is assumed that the bathymetry follows +! a normal distribution with sigma_b = 2.5d0. An improvement would +! be to provide the distribution based on high resolution data. +! +! Dupont, F. Dumont, D., Lemieux, J.F., Dumas-Lefebvre, E., Caya, A. +! in prep. +! +! authors: D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, F. Dupont +! + subroutine seabed_stress_factor_prob (nx_block, ny_block, & + icellt, indxti, indxtj, & + icellu, indxui, indxuj, & + aicen, vicen, & + hwater, Tbu) +! use modules + + use ice_arrays_column, only: hin_max + use ice_domain_size, only: ncat + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icellt, icellu ! no. of cells where icetmask = 1 + + integer (kind=int_kind), dimension (nx_block*ny_block), & + intent(in) :: & + indxti , & ! compressed index in i-direction + indxtj , & ! compressed index in j-direction + indxui , & ! compressed index in i-direction + indxuj ! compressed index in j-direction + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + hwater ! water depth at tracer location (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(in) :: & + aicen, & ! partial concentration for last thickness category in ITD + vicen ! partial volume for last thickness category in ITD (m) + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & + Tbu ! seabed stress factor (N/m^2) + +! local variables + + integer (kind=int_kind) :: & + i, j, ij, ii, n + + integer, parameter :: & + ncat_b = 100, & ! number of bathymetry categories + ncat_i = 100 ! number of ice thickness categories (log-normal) + + real (kind=dbl_kind), parameter :: & + max_depth = 50.0_dbl_kind, & ! initial range of log-normal distribution + mu_s = 0.1_dbl_kind, & ! friction coefficient + sigma_b = 2.5d0 ! Standard deviation of bathymetry + + real (kind=dbl_kind), dimension(ncat_i) :: & ! log-normal for ice thickness + x_k, & ! center of thickness categories (m) + g_k, & ! probability density function (thickness, 1/m) + P_x ! probability for each thickness category + + real (kind=dbl_kind), dimension(ncat_b) :: & ! normal dist for bathymetry + y_n, & ! center of bathymetry categories (m) + b_n, & ! probability density function (bathymetry, 1/m) + P_y ! probability for each bathymetry category + + real (kind=dbl_kind), dimension(ncat) :: & + vcat, acat + + integer, dimension(ncat_b) :: & + tmp ! Temporary vector tmp = merge(1,0,gt) + + logical, dimension (ncat_b) :: & + gt + + real (kind=dbl_kind) :: wid_i, wid_b, mu_i, sigma_i, mu_b, m_i, v_i ! parameters for PDFs + real (kind=dbl_kind), dimension(ncat_i):: tb_tmp + real (kind=dbl_kind), dimension (nx_block,ny_block):: Tbt ! seabed stress factor at t point (N/m^2) + real (kind=dbl_kind) :: atot, x_kmax + real (kind=dbl_kind) :: cut, rhoi, rhow, gravit, pi, puny + + character(len=*), parameter :: subname = '(seabed2_stress_coeff)' + call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi) + call icepack_query_parameters(gravit_out=gravit) + call icepack_query_parameters(pi_out=pi) + call icepack_query_parameters(puny_out=puny) + + Tbt=c0 + + do ij = 1, icellt + i = indxti(ij) + j = indxtj(ij) + + atot = sum(aicen(i,j,1:ncat)) + + if (atot > 0.05_dbl_kind .and. hwater(i,j) < max_depth) then + + mu_b = hwater(i,j) ! mean of PDF (normal dist) bathymetry + wid_i = max_depth/ncat_i ! width of ice categories + wid_b = c6*sigma_b/ncat_b ! width of bathymetry categories (6 sigma_b = 2x3 sigma_b) + + x_k = (/( wid_i*( real(i,kind=dbl_kind) - p5 ), i=1, ncat_i )/) + y_n = (/( ( mu_b-c3*sigma_b )+( real(i,kind=dbl_kind) - p5 )*( c6*sigma_b/ncat_b ), i=1, ncat_b )/) + + vcat(1:ncat) = vicen(i,j,1:ncat) + acat(1:ncat) = aicen(i,j,1:ncat) + + m_i = sum(vcat) + + v_i=c0 + do n =1, ncat + v_i = v_i + vcat(n)**2 / (max(acat(n), puny)) + enddo + v_i = v_i - m_i**2 + + mu_i = log(m_i/sqrt(c1 + v_i/m_i**2)) ! parameters for the log-normal + sigma_i = sqrt(log(c1 + v_i/m_i**2)) + + ! max thickness associated with percentile of log-normal PDF + ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al.) + + x_kmax = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0) + + ! Set x_kmax to hlev of the last category where there is ice + ! when there is no ice in the last category + cut = x_k(ncat_i) + do n = ncat,-1,1 + if (acat(n) < puny) then + cut = hin_max(n-1) + else + exit + endif + enddo + x_kmax = min(cut, x_kmax) + + g_k = exp(-(log(x_k) - mu_i) ** 2 / (c2 * sigma_i ** 2)) / (x_k * sigma_i * sqrt(c2 * pi)) + + b_n = exp(-(y_n - mu_b) ** 2 / (c2 * sigma_b ** 2)) / (sigma_b * sqrt(c2 * pi)) + + P_x = g_k*wid_i + P_y = b_n*wid_b + + do n =1, ncat_i + if (x_k(n) > x_kmax) P_x(n)=c0 + enddo + + ! calculate Tb factor at t-location + do n=1, ncat_i + gt = (y_n <= rhoi*x_k(n)/rhow) + tmp = merge(1,0,gt) + ii = sum(tmp) + if (ii == 0) then + tb_tmp(n) = c0 + else + tb_tmp(n) = max(mu_s*gravit*P_x(n)*sum(P_y(1:ii)*(rhoi*x_k(n) - rhow*y_n(1:ii))),c0) + endif + enddo + Tbt(i,j) = sum(tb_tmp)*exp(-alphab * (c1 - atot)) + endif + enddo + + do ij = 1, icellu + i = indxui(ij) + j = indxuj(ij) + ! convert quantities to u-location + Tbu(i,j) = max(Tbt(i,j),Tbt(i+1,j),Tbt(i,j+1),Tbt(i+1,j+1)) + enddo ! ij + + end subroutine seabed_stress_factor_prob + !======================================================================= ! Computes principal stresses for comparison with the theoretical diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 570e202c2..457a73ade 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -46,9 +46,9 @@ 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, basal_stress_coeff, basalstress, Ktens, & - stack_velocity_field, unstack_velocity_field + ecci, 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 use ice_flux, only: fm use ice_global_reductions, only: global_sum, global_allreduce_sum @@ -436,17 +436,31 @@ subroutine implicit_solver (dt) endif !----------------------------------------------------------------- - ! basal stress coefficients (landfast ice) + ! seabed stress factor Tbu (Tbu is part of Cb coefficient) !----------------------------------------------------------------- - if (basalstress) then + if (seabed_stress) then + !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks - call basal_stress_coeff (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:,iblk), & - vice (:,:,iblk), aice(:,:,iblk), & - hwater(:,:,iblk), Tbu (:,:,iblk)) + + if ( seabed_stress_method == 'LKD' ) then + + call seabed_stress_factor_LKD (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + + elseif ( seabed_stress_method == 'probabilistic' ) then + + call seabed_stress_factor_prob (nx_block, ny_block, & + icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & + icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + aicen(:,:,:,iblk), vicen(:,:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + endif + enddo !$OMP END PARALLEL DO endif @@ -527,7 +541,7 @@ subroutine implicit_solver (dt) !----------------------------------------------------------------- ! Compute seabed stress (diagnostic) !----------------------------------------------------------------- - if (basalstress) then + if (seabed_stress) then !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_seabed_stress (nx_block , ny_block , & @@ -1406,7 +1420,7 @@ subroutine calc_vrel_Cb (nx_block, ny_block, & uvel , vvel , & vrel , Cb) - use ice_dyn_shared, only: u0 ! residual velocity for basal stress (m/s) + use ice_dyn_shared, only: u0 ! residual velocity for seabed stress (m/s) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1417,7 +1431,7 @@ subroutine calc_vrel_Cb (nx_block, ny_block, & indxuj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & - Tbu, & ! coefficient for basal stress (N/m^2) + Tbu, & ! seabed stress factor (N/m^2) aiu , & ! ice fraction on u-grid uocn , & ! ocean current, x-direction (m/s) vocn , & ! ocean current, y-direction (m/s) @@ -1552,7 +1566,7 @@ subroutine matvec (nx_block, ny_block, & uvel , & ! x-component of velocity (m/s) vvel , & ! y-component of velocity (m/s) vrel , & ! coefficient for tauw - Cb , & ! coefficient for basal stress + Cb , & ! coefficient for seabed stress umassdti, & ! mass of U-cell/dt (kg/m^2 s) fm , & ! Coriolis param. * mass in U-cell (kg/s) uarear ! 1/uarea @@ -2361,7 +2375,7 @@ subroutine formDiag_step2 (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & vrel, & ! coefficient for tauw - Cb, & ! coefficient for basal stress + Cb, & ! coefficient for seabed stress umassdti, & ! mass of U-cell/dt (kg/m^2 s) uarear ! 1/uarea diff --git a/cicecore/cicedynB/dynamics/ice_transport_driver.F90 b/cicecore/cicedynB/dynamics/ice_transport_driver.F90 index 82e04dc71..c500e1631 100644 --- a/cicecore/cicedynB/dynamics/ice_transport_driver.F90 +++ b/cicecore/cicedynB/dynamics/ice_transport_driver.F90 @@ -1,7 +1,6 @@ !======================================================================= ! -!deprecate upwind Drivers for remapping and upwind ice transport -! Drivers for incremental remapping ice transport +! Drivers for remapping and upwind ice transport ! ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL ! @@ -10,7 +9,6 @@ ! 2006: Incorporated remap transport driver and renamed from ! ice_transport_upwind. ! 2011: ECH moved edgearea arrays into ice_transport_remap.F90 -! 2020: deprecated upwind transport module ice_transport_driver @@ -30,13 +28,12 @@ module ice_transport_driver implicit none private - public :: init_transport, transport_remap!deprecate upwind:, transport_upwind + public :: init_transport, transport_remap, transport_upwind character (len=char_len), public :: & advection ! type of advection scheme used -!deprecate upwind ! 'upwind' => 1st order donor cell scheme + ! 'upwind' => 1st order donor cell scheme ! 'remap' => remapping scheme - ! 'none' => advection off (ktransport = -1 also turns it off) logical, parameter :: & ! if true, prescribe area flux across each edge l_fixed_area = .false. @@ -72,9 +69,8 @@ module ice_transport_driver !======================================================================= ! ! This subroutine is a wrapper for init_remap, which initializes the -! remapping transport scheme. -!deprecate upwind If the model is run with upwind -!deprecate upwind! transport, no initializations are necessary. +! remapping transport scheme. If the model is run with upwind +! transport, no initializations are necessary. ! ! authors William H. Lipscomb, LANL @@ -684,12 +680,11 @@ subroutine transport_remap (dt) end subroutine transport_remap !======================================================================= -!deprecate upwind! +! ! Computes the transport equations for one timestep using upwind. Sets ! several fields into a work array and passes it to upwind routine. -!deprecate upwind - subroutine transport_upwind_deprecated (dt) + subroutine transport_upwind (dt) use ice_boundary, only: ice_HaloUpdate use ice_blocks, only: nx_block, ny_block, block, get_block, nx_block, ny_block @@ -774,52 +769,52 @@ subroutine transport_upwind_deprecated (dt) field_loc_Nface, field_type_vector) call ice_timer_stop(timer_bound) -!deprecate upwind !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) -!deprecate upwind do iblk = 1, nblocks -!deprecate upwind this_block = get_block(blocks_ice(iblk),iblk) -!deprecate upwind ilo = this_block%ilo -!deprecate upwind ihi = this_block%ihi -!deprecate upwind jlo = this_block%jlo -!deprecate upwind jhi = this_block%jhi + !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + !----------------------------------------------------------------- ! fill work arrays with fields to be advected !----------------------------------------------------------------- -!deprecate upwind -!deprecate upwind call state_to_work (nx_block, ny_block, & -!deprecate upwind ntrcr, & -!deprecate upwind narr, trcr_depend, & -!deprecate upwind aicen (:,:, :,iblk), trcrn (:,:,:,:,iblk), & -!deprecate upwind vicen (:,:, :,iblk), vsnon (:,:, :,iblk), & -!deprecate upwind aice0 (:,:, iblk), works (:,:, :,iblk)) + call state_to_work (nx_block, ny_block, & + ntrcr, & + narr, trcr_depend, & + aicen (:,:, :,iblk), trcrn (:,:,:,:,iblk), & + vicen (:,:, :,iblk), vsnon (:,:, :,iblk), & + aice0 (:,:, iblk), works (:,:, :,iblk)) !----------------------------------------------------------------- ! advect !----------------------------------------------------------------- -!deprecate upwind call upwind_field (nx_block, ny_block, & -!deprecate upwind ilo, ihi, jlo, jhi, & -!deprecate upwind dt, & -!deprecate upwind narr, works(:,:,:,iblk), & -!deprecate upwind uee(:,:,iblk), vnn (:,:,iblk), & -!deprecate upwind HTE(:,:,iblk), HTN (:,:,iblk), & -!deprecate upwind tarea(:,:,iblk)) + call upwind_field (nx_block, ny_block, & + ilo, ihi, jlo, jhi, & + dt, & + narr, works(:,:,:,iblk), & + uee(:,:,iblk), vnn (:,:,iblk), & + HTE(:,:,iblk), HTN (:,:,iblk), & + tarea(:,:,iblk)) !----------------------------------------------------------------- ! convert work arrays back to state variables !----------------------------------------------------------------- -!deprecate upwind call work_to_state (nx_block, ny_block, & -!deprecate upwind ntrcr, narr, & -!deprecate upwind trcr_depend(:), trcr_base(:,:), & -!deprecate upwind n_trcr_strata(:), nt_strata(:,:), & -!deprecate upwind aicen(:,:, :,iblk), trcrn (:,:,:,:,iblk), & -!deprecate upwind vicen(:,:, :,iblk), vsnon (:,:, :,iblk), & -!deprecate upwind aice0(:,:, iblk), works (:,:, :,iblk)) + call work_to_state (nx_block, ny_block, & + ntrcr, narr, & + trcr_depend(:), trcr_base(:,:), & + n_trcr_strata(:), nt_strata(:,:), & + aicen(:,:, :,iblk), trcrn (:,:,:,:,iblk), & + vicen(:,:, :,iblk), vsnon (:,:, :,iblk), & + aice0(:,:, iblk), works (:,:, :,iblk)) -!deprecate upwind enddo ! iblk -!deprecate upwind !$OMP END PARALLEL DO + enddo ! iblk + !$OMP END PARALLEL DO deallocate (works) @@ -837,8 +832,7 @@ subroutine transport_upwind_deprecated (dt) call ice_timer_stop(timer_advect) ! advection - end subroutine transport_upwind_deprecated -!deprecate upwind + end subroutine transport_upwind !======================================================================= ! The next few subroutines (through check_monotonicity) are called @@ -1461,12 +1455,12 @@ subroutine check_monotonicity (nx_block, ny_block, & end subroutine check_monotonicity !======================================================================= -!deprecate upwind! The remaining subroutines are called by transport_upwind. +! The remaining subroutines are called by transport_upwind. !======================================================================= ! ! Fill work array with state variables in preparation for upwind transport -!deprecate upwind - subroutine state_to_work_deprecated (nx_block, ny_block, & + + subroutine state_to_work (nx_block, ny_block, & ntrcr, & narr, trcr_depend, & aicen, trcrn, & @@ -1607,13 +1601,13 @@ subroutine state_to_work_deprecated (nx_block, ny_block, & if (narr /= narrays) write(nu_diag,*) & "Wrong number of arrays in transport bound call" - end subroutine state_to_work_deprecated + end subroutine state_to_work !======================================================================= ! ! Convert work array back to state variables -!deprecate upwind - subroutine work_to_state_deprecated (nx_block, ny_block, & + + subroutine work_to_state (nx_block, ny_block, & ntrcr, narr, & trcr_depend, & trcr_base, & @@ -1721,13 +1715,13 @@ subroutine work_to_state_deprecated (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - end subroutine work_to_state_deprecated + end subroutine work_to_state !======================================================================= ! ! upwind transport algorithm -!deprecate upwind - subroutine upwind_field_deprecated (nx_block, ny_block, & + + subroutine upwind_field (nx_block, ny_block, & ilo, ihi, jlo, jhi, & dt, & narrays, phi, & @@ -1770,26 +1764,26 @@ subroutine upwind_field_deprecated (nx_block, ny_block, & do n = 1, narrays -!deprecate upwind do j = 1, jhi -!deprecate upwind do i = 1, ihi -!deprecate upwind worka(i,j)= & -!deprecate upwind upwind(phi(i,j,n),phi(i+1,j,n),uee(i,j),HTE(i,j),dt) -!deprecate upwind workb(i,j)= & -!deprecate upwind upwind(phi(i,j,n),phi(i,j+1,n),vnn(i,j),HTN(i,j),dt) -!deprecate upwind enddo -!deprecate upwind enddo - -!deprecate upwind do j = jlo, jhi -!deprecate upwind do i = ilo, ihi -!deprecate upwind phi(i,j,n) = phi(i,j,n) - ( worka(i,j)-worka(i-1,j) & -!deprecate upwind + workb(i,j)-workb(i,j-1) ) & -!deprecate upwind / tarea(i,j) -!deprecate upwind enddo -!deprecate upwind enddo + do j = 1, jhi + do i = 1, ihi + worka(i,j)= & + upwind(phi(i,j,n),phi(i+1,j,n),uee(i,j),HTE(i,j),dt) + workb(i,j)= & + upwind(phi(i,j,n),phi(i,j+1,n),vnn(i,j),HTN(i,j),dt) + enddo + enddo + + do j = jlo, jhi + do i = ilo, ihi + phi(i,j,n) = phi(i,j,n) - ( worka(i,j)-worka(i-1,j) & + + workb(i,j)-workb(i,j-1) ) & + / tarea(i,j) + enddo + enddo enddo ! narrays - end subroutine upwind_field_deprecated + end subroutine upwind_field !======================================================================= @@ -1797,13 +1791,13 @@ end subroutine upwind_field_deprecated ! Define upwind function !------------------------------------------------------------------- -!deprecate upwind real(kind=dbl_kind) function upwind(y1,y2,a,h,dt) + real(kind=dbl_kind) function upwind(y1,y2,a,h,dt) -!deprecate upwind real(kind=dbl_kind), intent(in) :: y1,y2,a,h,dt + real(kind=dbl_kind), intent(in) :: y1,y2,a,h,dt -!deprecate upwind upwind = p5*dt*h*((a+abs(a))*y1+(a-abs(a))*y2) + upwind = p5*dt*h*((a+abs(a))*y1+(a-abs(a))*y2) -!deprecate upwind end function upwind + end function upwind !======================================================================= diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index 97b726fdb..71253a4b1 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -46,7 +46,7 @@ module ice_flux vocn , & ! ocean current, y-direction (m/s) ss_tltx , & ! sea surface slope, x-direction (m/m) ss_tlty , & ! sea surface slope, y-direction - hwater , & ! water depth for basal stress calc (landfast ice) + hwater , & ! water depth for seabed stress calc (landfast ice) ! out to atmosphere strairxT, & ! stress on ice by air, x-direction @@ -63,8 +63,8 @@ module ice_flux sig1 , & ! normalized principal stress component sig2 , & ! normalized principal stress component sigP , & ! internal ice pressure (N/m) - taubx , & ! basal stress (x) (N/m^2) - tauby , & ! basal stress (y) (N/m^2) + taubx , & ! seabed stress (x) (N/m^2) + tauby , & ! seabed stress (y) (N/m^2) strairx , & ! stress on ice by air, x-direction strairy , & ! stress on ice by air, y-direction strocnx , & ! ice-ocean stress, x-direction @@ -112,7 +112,7 @@ module ice_flux real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & fm , & ! Coriolis param. * mass in U-cell (kg/s) - Tbu ! coefficient for basal stress (N/m^2) + Tbu ! factor for seabed stress (N/m^2) !----------------------------------------------------------------- ! Thermodynamic component @@ -351,7 +351,7 @@ subroutine alloc_flux vocn (nx_block,ny_block,max_blocks), & ! ocean current, y-direction (m/s) ss_tltx (nx_block,ny_block,max_blocks), & ! sea surface slope, x-direction (m/m) ss_tlty (nx_block,ny_block,max_blocks), & ! sea surface slope, y-direction - hwater (nx_block,ny_block,max_blocks), & ! water depth for basal stress calc (landfast ice) + hwater (nx_block,ny_block,max_blocks), & ! water depth for seabed stress calc (landfast ice) strairxT (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction strairyT (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction strocnxT (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction @@ -359,8 +359,8 @@ subroutine alloc_flux sig1 (nx_block,ny_block,max_blocks), & ! normalized principal stress component sig2 (nx_block,ny_block,max_blocks), & ! normalized principal stress component sigP (nx_block,ny_block,max_blocks), & ! internal ice pressure (N/m) - taubx (nx_block,ny_block,max_blocks), & ! basal stress (x) (N/m^2) - tauby (nx_block,ny_block,max_blocks), & ! basal stress (y) (N/m^2) + taubx (nx_block,ny_block,max_blocks), & ! seabed stress (x) (N/m^2) + tauby (nx_block,ny_block,max_blocks), & ! seabed stress (y) (N/m^2) strairx (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction strairy (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction strocnx (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction @@ -390,7 +390,7 @@ subroutine alloc_flux stress12_4 (nx_block,ny_block,max_blocks), & ! sigma12 iceumask (nx_block,ny_block,max_blocks), & ! ice extent mask (U-cell) fm (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in U-cell (kg/s) - Tbu (nx_block,ny_block,max_blocks), & ! coefficient for basal stress (landfast ice) + Tbu (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice) zlvl (nx_block,ny_block,max_blocks), & ! atm level height (m) uatm (nx_block,ny_block,max_blocks), & ! wind velocity components (m/s) vatm (nx_block,ny_block,max_blocks), & diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 old mode 100755 new mode 100644 index 43cf92a48..edbba8101 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -533,17 +533,21 @@ subroutine get_forcing_atmo integer (kind=int_kind) :: & iblk, & ! block index ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + fyear_old, & ! prior fyear value nt_Tsfc type (block) :: & this_block ! block information for current block character(len=*), parameter :: subname = '(get_forcing_atmo)' - + + fyear_old = fyear fyear = fyear_init + mod(nyr-1,ycycle) ! current year - if (trim(atm_data_type) /= 'default' .and. istep <= 1 & - .and. my_task == master_task) then - write (nu_diag,*) ' Current forcing data year = ',fyear + if (trim(atm_data_type) /= 'default' .and. & + (istep <= 1 .or. fyear /= fyear_old)) then + if (my_task == master_task) then + write (nu_diag,*) ' Current forcing data year = ',fyear + endif endif call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc) @@ -2322,162 +2326,174 @@ subroutine JRA55_data (yr) use ice_state, only: aice use ice_calendar, only: days_per_year, use_leap_years + integer (kind=int_kind), intent(in) :: & + yr ! current forcing year + integer (kind=int_kind) :: & ncid , & ! netcdf file id - i, j , & - ixm,ixx,ixp , & ! record numbers for neighboring months + i, j, n1, iblk, & + yrp , & ! year after yr in forcing cycle recnum , & ! record number maxrec , & ! maximum record number recslot , & ! spline slot for current record - midmonth , & ! middle day of month - dataloc , & ! = 1 for data located in middle of time interval + dataloc ! = 1 for data located in middle of time interval ! = 2 for date located at end of time interval - iblk , & ! block index - ilo,ihi,jlo,jhi, & ! beginning and end of physical domain - yr ! current forcing year real (kind=dbl_kind) :: & sec3hr , & ! number of seconds in 3 hours secday , & ! number of seconds in day + eps, tt , & ! interpolation coeff calc Tffresh , & vmin, vmax - logical (kind=log_kind) :: readm, read6,debug_n_d - - type (block) :: & - this_block ! block information for current block + logical (kind=log_kind) :: debug_n_d = .false. + character (char_len_long) :: uwind_file_old character(len=64) :: fieldname !netcdf field name character(len=*), parameter :: subname = '(JRA55_data)' - debug_n_d = .false. !usually false - call icepack_query_parameters(Tffresh_out=Tffresh) call icepack_query_parameters(secday_out=secday) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - !------------------------------------------------------------------- - ! 3-hourly data - ! - ! Assume that the 3-hourly value is located at the end of the - ! 3-hour period. This is the convention for NCEP reanalysis data. - ! E.g. record 1 gives conditions at 3 am GMT on 1 January. - !------------------------------------------------------------------- - - dataloc = 2 ! data located at end of interval sec3hr = secday/c8 ! seconds in 3 hours - !maxrec = 2920 ! 365*8; for leap years = 366*8 - - if (use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar maxrec = days_per_year*8 - if(days_per_year == 365 .and. (mod(yr, 4) == 0)) then - call abort_ice('days_per_year should be set to 366 for leap years') - end if + if (debug_n_d .and. my_task == master_task) then + write (nu_diag,*) subname,'recnum',recnum + write (nu_diag,*) subname,'maxrec',maxrec + write (nu_diag,*) subname,'days_per_year', days_per_year + endif - ! current record number - recnum = 8*int(yday) - 7 + int(real(sec,kind=dbl_kind)/sec3hr) + !------------------------------------------------------------------- + ! 3-hourly data + ! states are instantaneous, 1st record is 00z Jan 1 + ! fluxes are 3 hour averages, 1st record is 00z-03z Jan 1 + ! Both states and fluxes have 1st record defined as 00z Jan 1 + ! interpolate states, do not interpolate fluxes + ! fluxes are held constant from [init period, end period) + !------------------------------------------------------------------- + ! File is NETCDF with winds in NORTH and EAST direction + ! file variable names are: + ! glbrad (shortwave W/m^2) + ! dlwsfc (longwave W/m^2) + ! wndewd (eastward wind m/s) + ! wndnwd (northward wind m/s) + ! airtmp (air temperature K) + ! spchmd (specific humidity kg/kg) + ! ttlpcp (precipitation kg/m s-1) + !------------------------------------------------------------------- - ! Compute record numbers for surrounding data (2 on each side) + uwind_file_old = uwind_file + call file_year(uwind_file,yr) + if (uwind_file /= uwind_file_old .and. my_task == master_task) then + write(nu_diag,*) subname,' reading forcing file = ',trim(uwind_file) + endif - ixm = mod(recnum+maxrec-2,maxrec) + 1 - ixx = mod(recnum-1, maxrec) + 1 + call ice_open_nc(uwind_file,ncid) - ! Compute interpolation coefficients - ! If data is located at the end of the time interval, then the - ! data value for the current record goes in slot 2 + do n1 = 1,2 - recslot = 2 - ixp = -99 - call interp_coeff (recnum, recslot, sec3hr, dataloc) + if (n1 == 1) then + recnum = 8*int(yday) - 7 + int(real(sec,kind=dbl_kind)/sec3hr) + if (my_task == master_task .and. (recnum <= 2 .or. recnum >= maxrec-1)) then + write(nu_diag,*) subname,' reading forcing file 1st ts = ',trim(uwind_file) + endif + elseif (n1 == 2) then + recnum = 8*int(yday) - 7 + int(real(sec,kind=dbl_kind)/sec3hr) + 1 + if (recnum > maxrec) then + yrp = fyear_init + mod(nyr,ycycle) ! next year + recnum = 1 + call file_year(uwind_file,yrp) + if (my_task == master_task) then + write(nu_diag,*) subname,' reading forcing file 2nd ts = ',trim(uwind_file) + endif + call ice_close_nc(ncid) + call ice_open_nc(uwind_file,ncid) + endif + endif - ! Read - read6 = .false. - if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true. - !------------------------------------------------------------------- - ! File is NETCDF with winds in NORTH and EAST direction - ! file variable names are: - ! glbrad (shortwave W/m^2) - ! dlwsfc (longwave W/m^2) - ! wndewd (eastward wind m/s) - ! wndnwd (northward wind m/s) - ! airtmp (air temperature K) - ! spchmd (specific humidity kg/kg) - ! ttlpcp (precipitation kg/m s-1) - !------------------------------------------------------------------- - call ice_open_nc(uwind_file,ncid) - - fieldname = 'airtmp' - call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,2,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) + if (debug_n_d .and. my_task == master_task) then + write(nu_diag,*) subname,' read recnum = ',recnum,n1 + endif - fieldname = 'wndewd' - call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,2,:),debug_n_d, & + fieldname = 'airtmp' + call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - fieldname = 'wndnwd' - call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,2,:),debug_n_d, & + fieldname = 'wndewd' + call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - fieldname = 'spchmd' - call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,2,:),debug_n_d, & + fieldname = 'wndnwd' + call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - fieldname = 'glbrad' - call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,2,:),debug_n_d, & + fieldname = 'spchmd' + call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - fieldname = 'dlwsfc' - call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,1,:),debug_n_d, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,2,:),debug_n_d, & + ! only read one timestep for fluxes, 3 hr average, no interpolation + if (n1 == 1) then + fieldname = 'glbrad' + call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - fieldname = 'ttlpcp' - call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,1,:),debug_n_d, & + fieldname = 'dlwsfc' + call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) - call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,2,:),debug_n_d, & + + fieldname = 'ttlpcp' + call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,n1,:),debug_n_d, & field_loc=field_loc_center, & field_type=field_type_scalar) + endif + + enddo + + call ice_close_nc(ncid) - call ice_close_nc(ncid) + ! reset uwind_file to original year + call file_year(uwind_file,yr) + ! Compute interpolation coefficients + eps = 1.0e-6 + tt = real(mod(sec,nint(sec3hr)),kind=dbl_kind) + c2intp = tt / sec3hr + if (c2intp < c0 .and. c2intp > c0-eps) c2intp = c0 + if (c2intp > c1 .and. c2intp < c1+eps) c2intp = c1 + c1intp = 1.0_dbl_kind - c2intp + if (c2intp < c0 .or. c2intp > c1) then + write(nu_diag,*) subname,' ERROR: c2intp = ',c2intp + call abort_ice (error_message=subname//' ERROR: c2intp out of range', & + file=__FILE__, line=__LINE__) + endif + if (debug_n_d .and. my_task == master_task) then + write(nu_diag,*) subname,' c12intp = ',c1intp,c2intp + endif ! Interpolate call interpolate_data (Tair_data, Tair) call interpolate_data (uatm_data, uatm) call interpolate_data (vatm_data, vatm) call interpolate_data (Qa_data, Qa) - call interpolate_data (fsw_data, fsw) - call interpolate_data (flw_data, flw) - call interpolate_data (fsnow_data, fsnow) + ! use 3 hr average for heat flux and precip fields + ! call interpolate_data (fsw_data, fsw) + ! call interpolate_data (flw_data, flw) + ! call interpolate_data (fsnow_data, fsnow) + fsw(:,:,:) = fsw_data(:,:,1,:) + flw(:,:,:) = flw_data(:,:,1,:) + fsnow(:,:,:) = fsnow_data(:,:,1,:) - !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks ! limit summer Tair values where ice is present do j = 1, ny_block @@ -2501,43 +2517,37 @@ subroutine JRA55_data (yr) enddo ! iblk !$OMP END PARALLEL DO - ! Save record number - oldrecnum = recnum - - if (dbug) then - if (my_task == master_task) write (nu_diag,*) 'JRA55_bulk_data' + if (debug_n_d .or. dbug) then + if (my_task.eq.master_task) & + write (nu_diag,*) subname,'JRA55_bulk_data' vmin = global_minval(fsw,distrb_info,tmask) - vmax = global_maxval(fsw,distrb_info,tmask) - if (my_task.eq.master_task) & - write (nu_diag,*) 'fsw',vmin,vmax + if (my_task.eq.master_task) & + write (nu_diag,*) subname,'fsw',vmin,vmax vmin = global_minval(flw,distrb_info,tmask) vmax = global_maxval(flw,distrb_info,tmask) if (my_task.eq.master_task) & - write (nu_diag,*) 'flw',vmin,vmax + write (nu_diag,*) subname,'flw',vmin,vmax vmin =global_minval(fsnow,distrb_info,tmask) vmax =global_maxval(fsnow,distrb_info,tmask) if (my_task.eq.master_task) & - write (nu_diag,*) 'fsnow',vmin,vmax + write (nu_diag,*) subname,'fsnow',vmin,vmax vmin = global_minval(Tair,distrb_info,tmask) vmax = global_maxval(Tair,distrb_info,tmask) if (my_task.eq.master_task) & - write (nu_diag,*) 'Tair',vmin,vmax + write (nu_diag,*) subname,'Tair',vmin,vmax vmin = global_minval(uatm,distrb_info,umask) vmax = global_maxval(uatm,distrb_info,umask) if (my_task.eq.master_task) & - write (nu_diag,*) 'uatm',vmin,vmax + write (nu_diag,*) subname,'uatm',vmin,vmax vmin = global_minval(vatm,distrb_info,umask) vmax = global_maxval(vatm,distrb_info,umask) if (my_task.eq.master_task) & - write (nu_diag,*) 'vatm',vmin,vmax + write (nu_diag,*) subname,'vatm',vmin,vmax vmin = global_minval(Qa,distrb_info,tmask) vmax = global_maxval(Qa,distrb_info,tmask) - if (my_task.eq.master_task) & - write (nu_diag,*) 'Qa',vmin,vmax - if (my_task.eq.master_task) & - write (nu_diag,*) 'maxrec',maxrec - write (nu_diag,*) 'days_per_year', days_per_year + if (my_task.eq.master_task) & + write (nu_diag,*) subname,'Qa',vmin,vmax endif ! dbug diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index fb9c45978..b59a93862 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -59,7 +59,7 @@ subroutine input_data use ice_broadcast, only: broadcast_scalar, broadcast_array use ice_diagnostics, only: diag_file, print_global, print_points, latpnt, lonpnt - use ice_domain, only: close_boundaries, ns_boundary_type, orca_halogrid + use ice_domain, only: close_boundaries, orca_halogrid use ice_domain_size, only: ncat, nilyr, nslyr, nblyr, nfsd, nfreq, & n_iso, n_aero, n_zaero, n_algae, & n_doc, n_dic, n_don, n_fed, n_fep, & @@ -97,9 +97,10 @@ subroutine input_data dxrect, dyrect use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, & kevp_kernel, & - basalstress, k1, k2, alphab, threshold_hw, & + seabed_stress, seabed_stress_method, & + k1, k2, alphab, threshold_hw, & Ktens, e_ratio, coriolis, ssh_stress, & - kridge, ktransport, brlx, arlx + 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, & @@ -124,10 +125,10 @@ subroutine input_data mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, & a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, & phi_c_slow_mode, phi_i_mushy, kalg, atmiter_conv, Pstar, Cstar, & - sw_frac, sw_dtemp + sw_frac, sw_dtemp, floediam, hfrazilmin integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, natmiter, & - kitd, kcatbound + kitd, kcatbound, ktransport character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & tfrz_option, frzpnd, atmbndy, wave_spec_type @@ -190,7 +191,7 @@ subroutine input_data kitd, ktherm, conduct, ksno, & a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, & dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy, & - sw_redist, sw_frac, sw_dtemp + floediam, hfrazilmin namelist /dynamics_nml/ & kdyn, ndte, revised_evp, yield_curve, & @@ -198,7 +199,7 @@ subroutine input_data brlx, arlx, ssh_stress, & advection, coriolis, kridge, ktransport, & kstrength, krdg_partic, krdg_redist, mu_rdg, & - e_ratio, Ktens, Cf, basalstress, & + e_ratio, Ktens, Cf, seabed_stress, & k1, maxits_nonlin, precond, dim_fgmres, & dim_pgmres, maxits_fgmres, maxits_pgmres, monitor_nonlin, & monitor_fgmres, monitor_pgmres, reltol_nonlin, reltol_fgmres, & @@ -206,12 +207,13 @@ subroutine input_data damping_andacc, start_andacc, fpfunc_andacc, use_mean_vrel, & ortho_type, & k2, alphab, threshold_hw, & - Pstar, Cstar + seabed_stress_method, Pstar, Cstar namelist /shortwave_nml/ & shortwave, albedo_type, & albicev, albicei, albsnowv, albsnowi, & ahmax, R_ice, R_pnd, R_snw, & + sw_redist, sw_frac, sw_dtemp, & dT_mlt, rsnw_mlt, kalg namelist /ponds_nml/ & @@ -308,7 +310,7 @@ subroutine input_data kitd = 1 ! type of itd conversions (0 = delta, 1 = linear) kcatbound = 1 ! category boundary formula (0 = old, 1 = new, etc) - kdyn = 1 ! type of dynamics (-1, 0 = off, 1 = evp, 2 = eap) + kdyn = 1 ! type of dynamics (-1, 0 = off, 1 = evp, 2 = eap, 3 = vp) ndtd = 1 ! dynamic time steps per thermodynamic time step ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte kevp_kernel = 0 ! EVP kernel (0 = 2D, >0: 1D. Only ver. 2 is implemented yet) @@ -327,9 +329,10 @@ subroutine input_data dxrect = 0.0_dbl_kind ! user defined grid spacing in cm in x direction dyrect = 0.0_dbl_kind ! user defined grid spacing in cm in y direction close_boundaries = .false. ! true = set land on edges of grid - basalstress= .false. ! if true, basal stress for landfast is on - k1 = 8.0_dbl_kind ! 1st free parameter for landfast parameterization - k2 = 15.0_dbl_kind ! dah: second free parameter (N/m^3) for landfast parametrization + seabed_stress= .false. ! if true, seabed stress for landfast is on + seabed_stress_method = 'LKD' ! LKD = Lemieux et al 2015, probabilistic = Dupont et al. in prep + k1 = 7.5_dbl_kind ! 1st free parameter for landfast parameterization + k2 = 15.0_dbl_kind ! 2nd free parameter (N/m^3) for landfast parametrization 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) @@ -367,7 +370,7 @@ subroutine input_data calc_Tsfc = .true. ! calculate surface temperature update_ocn_f = .false. ! include fresh water and salt fluxes for frazil ustar_min = 0.005 ! minimum friction velocity for ocean heat flux (m/s) - emissivity = 0.95 ! emissivity of snow and ice + emissivity = 0.985 ! emissivity of snow and ice l_mpond_fresh = .false. ! logical switch for including meltpond freshwater ! flux feedback to ocean model fbot_xfer_type = 'constant' ! transfer coefficient type for ocn heat flux @@ -473,6 +476,9 @@ subroutine input_data phi_c_slow_mode = 0.05_dbl_kind ! critical liquid fraction porosity cutoff phi_i_mushy = 0.85_dbl_kind ! liquid fraction of congelation ice + floediam = 300.0_dbl_kind ! min thickness of new frazil ice (m) + hfrazilmin = 0.05_dbl_kind ! effective floe diameter (m) + ! shortwave redistribution in the thermodynamics sw_redist = .false. sw_frac = 0.9_dbl_kind @@ -573,218 +579,221 @@ subroutine input_data ! broadcast namelist settings !----------------------------------------------------------------- - call broadcast_scalar(numin, master_task) - call broadcast_scalar(numax, master_task) - call broadcast_scalar(days_per_year, master_task) - call broadcast_scalar(use_leap_years, master_task) - call broadcast_scalar(year_init, master_task) - call broadcast_scalar(istep0, master_task) - call broadcast_scalar(dt, master_task) - call broadcast_scalar(npt, master_task) - call broadcast_scalar(diagfreq, master_task) - call broadcast_scalar(print_points, master_task) - call broadcast_scalar(print_global, master_task) - call broadcast_scalar(bfbflag, master_task) - call broadcast_scalar(diag_type, master_task) - call broadcast_scalar(diag_file, master_task) + call broadcast_scalar(numin, master_task) + call broadcast_scalar(numax, master_task) + call broadcast_scalar(days_per_year, master_task) + call broadcast_scalar(use_leap_years, master_task) + call broadcast_scalar(year_init, master_task) + call broadcast_scalar(istep0, master_task) + call broadcast_scalar(dt, master_task) + call broadcast_scalar(npt, master_task) + call broadcast_scalar(diagfreq, master_task) + call broadcast_scalar(print_points, master_task) + call broadcast_scalar(print_global, master_task) + call broadcast_scalar(bfbflag, master_task) + call broadcast_scalar(diag_type, master_task) + call broadcast_scalar(diag_file, master_task) do n = 1, max_nstrm - call broadcast_scalar(histfreq(n), master_task) + call broadcast_scalar(histfreq(n), master_task) enddo - call broadcast_array(histfreq_n, master_task) - call broadcast_scalar(hist_avg, master_task) - call broadcast_scalar(history_dir, master_task) - call broadcast_scalar(history_file, master_task) - call broadcast_scalar(history_precision, master_task) - call broadcast_scalar(history_format, master_task) - call broadcast_scalar(write_ic, master_task) - call broadcast_scalar(cpl_bgc, master_task) - call broadcast_scalar(incond_dir, master_task) - call broadcast_scalar(incond_file, master_task) - call broadcast_scalar(dumpfreq, master_task) - call broadcast_scalar(dumpfreq_n, master_task) - call broadcast_scalar(dump_last, master_task) - call broadcast_scalar(restart_file, master_task) - call broadcast_scalar(restart, master_task) - call broadcast_scalar(restart_dir, master_task) - call broadcast_scalar(restart_ext, master_task) - call broadcast_scalar(restart_coszen, master_task) - call broadcast_scalar(use_restart_time, master_task) - call broadcast_scalar(restart_format, master_task) - call broadcast_scalar(lcdf64, master_task) - call broadcast_scalar(pointer_file, master_task) - call broadcast_scalar(ice_ic, master_task) - call broadcast_scalar(grid_format, master_task) - call broadcast_scalar(dxrect, master_task) - call broadcast_scalar(dyrect, master_task) - call broadcast_scalar(close_boundaries, master_task) - call broadcast_scalar(grid_type, master_task) - call broadcast_scalar(grid_file, master_task) - call broadcast_scalar(gridcpl_file, master_task) - call broadcast_scalar(orca_halogrid, master_task) - call broadcast_scalar(bathymetry_file, master_task) - call broadcast_scalar(bathymetry_format, master_task) - call broadcast_scalar(use_bathymetry, master_task) - call broadcast_scalar(kmt_file, master_task) - call broadcast_scalar(kitd, master_task) - call broadcast_scalar(kcatbound, master_task) - call broadcast_scalar(kdyn, master_task) - call broadcast_scalar(ndtd, master_task) - call broadcast_scalar(ndte, master_task) - call broadcast_scalar(kevp_kernel, master_task) - call broadcast_scalar(brlx, master_task) - call broadcast_scalar(arlx, master_task) - call broadcast_scalar(revised_evp, master_task) - call broadcast_scalar(yield_curve, master_task) - call broadcast_scalar(kstrength, master_task) - call broadcast_scalar(Pstar, master_task) - call broadcast_scalar(Cstar, master_task) - call broadcast_scalar(krdg_partic, master_task) - call broadcast_scalar(krdg_redist, master_task) - call broadcast_scalar(mu_rdg, master_task) - call broadcast_scalar(Cf, master_task) - call broadcast_scalar(ksno, master_task) - call broadcast_scalar(basalstress, master_task) - call broadcast_scalar(k1, master_task) - call broadcast_scalar(k2, master_task) - 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(advection, master_task) - call broadcast_scalar(conserv_check, master_task) - call broadcast_scalar(shortwave, master_task) - call broadcast_scalar(albedo_type, master_task) - call broadcast_scalar(ktherm, master_task) - call broadcast_scalar(coriolis, master_task) - call broadcast_scalar(ssh_stress, master_task) - call broadcast_scalar(kridge, master_task) - call broadcast_scalar(ktransport, master_task) - call broadcast_scalar(maxits_nonlin, master_task) - call broadcast_scalar(precond, master_task) - call broadcast_scalar(dim_fgmres, master_task) - call broadcast_scalar(dim_pgmres, master_task) - call broadcast_scalar(maxits_fgmres, master_task) - call broadcast_scalar(maxits_pgmres, master_task) - call broadcast_scalar(monitor_nonlin, master_task) - call broadcast_scalar(monitor_fgmres, master_task) - call broadcast_scalar(monitor_pgmres, master_task) - call broadcast_scalar(ortho_type, master_task) - call broadcast_scalar(reltol_nonlin, master_task) - call broadcast_scalar(reltol_fgmres, master_task) - call broadcast_scalar(reltol_pgmres, master_task) - call broadcast_scalar(algo_nonlin, master_task) - call broadcast_scalar(fpfunc_andacc, master_task) - call broadcast_scalar(dim_andacc, master_task) - call broadcast_scalar(reltol_andacc, master_task) - call broadcast_scalar(damping_andacc, master_task) - call broadcast_scalar(start_andacc, master_task) - call broadcast_scalar(use_mean_vrel, master_task) - call broadcast_scalar(conduct, master_task) - call broadcast_scalar(R_ice, master_task) - call broadcast_scalar(R_pnd, master_task) - call broadcast_scalar(R_snw, master_task) - call broadcast_scalar(dT_mlt, master_task) - call broadcast_scalar(rsnw_mlt, master_task) - call broadcast_scalar(kalg, master_task) - call broadcast_scalar(hp1, master_task) - call broadcast_scalar(hs0, master_task) - call broadcast_scalar(hs1, master_task) - call broadcast_scalar(dpscale, master_task) - call broadcast_scalar(frzpnd, master_task) - call broadcast_scalar(rfracmin, master_task) - call broadcast_scalar(rfracmax, master_task) - call broadcast_scalar(pndaspect, master_task) - call broadcast_scalar(albicev, master_task) - call broadcast_scalar(albicei, master_task) - call broadcast_scalar(albsnowv, master_task) - call broadcast_scalar(albsnowi, master_task) - call broadcast_scalar(ahmax, master_task) - call broadcast_scalar(atmbndy, master_task) - call broadcast_scalar(fyear_init, master_task) - call broadcast_scalar(ycycle, master_task) - call broadcast_scalar(atm_data_format, master_task) - call broadcast_scalar(atm_data_type, master_task) - call broadcast_scalar(atm_data_dir, master_task) - call broadcast_scalar(rotate_wind, master_task) - call broadcast_scalar(calc_strair, master_task) - call broadcast_scalar(calc_Tsfc, master_task) - call broadcast_scalar(formdrag, master_task) - call broadcast_scalar(highfreq, master_task) - call broadcast_scalar(natmiter, master_task) - call broadcast_scalar(atmiter_conv, master_task) - call broadcast_scalar(update_ocn_f, master_task) - call broadcast_scalar(l_mpond_fresh, master_task) - call broadcast_scalar(ustar_min, master_task) - call broadcast_scalar(emissivity, master_task) - call broadcast_scalar(fbot_xfer_type, master_task) - call broadcast_scalar(precip_units, master_task) - call broadcast_scalar(oceanmixed_ice, master_task) - call broadcast_scalar(wave_spec_type, master_task) - call broadcast_scalar(wave_spec_file, master_task) - call broadcast_scalar(nfreq, master_task) - call broadcast_scalar(tfrz_option, master_task) - call broadcast_scalar(ocn_data_format, master_task) - call broadcast_scalar(bgc_data_type, master_task) - call broadcast_scalar(fe_data_type, master_task) - call broadcast_scalar(ice_data_type, master_task) - call broadcast_scalar(bgc_data_dir, master_task) - call broadcast_scalar(ocn_data_type, master_task) - call broadcast_scalar(ocn_data_dir, master_task) - call broadcast_scalar(oceanmixed_file, master_task) - call broadcast_scalar(restore_ocn, master_task) - call broadcast_scalar(trestore, master_task) - call broadcast_scalar(restore_ice, master_task) - call broadcast_scalar(dbug, master_task) - call broadcast_array (latpnt(1:2), master_task) - call broadcast_array (lonpnt(1:2), master_task) - call broadcast_scalar(runid, master_task) - call broadcast_scalar(runtype, master_task) + call broadcast_array(histfreq_n, master_task) + call broadcast_scalar(hist_avg, master_task) + call broadcast_scalar(history_dir, master_task) + call broadcast_scalar(history_file, master_task) + call broadcast_scalar(history_precision, master_task) + call broadcast_scalar(history_format, master_task) + call broadcast_scalar(write_ic, master_task) + call broadcast_scalar(cpl_bgc, master_task) + call broadcast_scalar(incond_dir, master_task) + call broadcast_scalar(incond_file, master_task) + call broadcast_scalar(dumpfreq, master_task) + call broadcast_scalar(dumpfreq_n, master_task) + call broadcast_scalar(dump_last, master_task) + call broadcast_scalar(restart_file, master_task) + call broadcast_scalar(restart, master_task) + call broadcast_scalar(restart_dir, master_task) + call broadcast_scalar(restart_ext, master_task) + call broadcast_scalar(restart_coszen, master_task) + call broadcast_scalar(use_restart_time, master_task) + call broadcast_scalar(restart_format, master_task) + call broadcast_scalar(lcdf64, master_task) + call broadcast_scalar(pointer_file, master_task) + call broadcast_scalar(ice_ic, master_task) + call broadcast_scalar(grid_format, master_task) + call broadcast_scalar(dxrect, master_task) + call broadcast_scalar(dyrect, master_task) + call broadcast_scalar(close_boundaries, master_task) + call broadcast_scalar(grid_type, master_task) + call broadcast_scalar(grid_file, master_task) + call broadcast_scalar(gridcpl_file, master_task) + call broadcast_scalar(orca_halogrid, master_task) + call broadcast_scalar(bathymetry_file, master_task) + call broadcast_scalar(bathymetry_format, master_task) + call broadcast_scalar(use_bathymetry, master_task) + call broadcast_scalar(kmt_file, master_task) + call broadcast_scalar(kitd, master_task) + call broadcast_scalar(kcatbound, master_task) + call broadcast_scalar(kdyn, master_task) + call broadcast_scalar(ndtd, master_task) + call broadcast_scalar(ndte, master_task) + call broadcast_scalar(kevp_kernel, master_task) + call broadcast_scalar(brlx, master_task) + call broadcast_scalar(arlx, master_task) + call broadcast_scalar(revised_evp, master_task) + call broadcast_scalar(yield_curve, master_task) + call broadcast_scalar(kstrength, master_task) + call broadcast_scalar(Pstar, master_task) + call broadcast_scalar(Cstar, master_task) + call broadcast_scalar(krdg_partic, master_task) + call broadcast_scalar(krdg_redist, master_task) + call broadcast_scalar(mu_rdg, master_task) + call broadcast_scalar(Cf, master_task) + call broadcast_scalar(ksno, master_task) + call broadcast_scalar(seabed_stress, master_task) + call broadcast_scalar(seabed_stress_method, master_task) + call broadcast_scalar(k1, master_task) + call broadcast_scalar(k2, master_task) + 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(advection, master_task) + call broadcast_scalar(conserv_check, master_task) + call broadcast_scalar(shortwave, master_task) + call broadcast_scalar(albedo_type, master_task) + call broadcast_scalar(ktherm, master_task) + call broadcast_scalar(coriolis, master_task) + call broadcast_scalar(ssh_stress, master_task) + call broadcast_scalar(kridge, master_task) + call broadcast_scalar(ktransport, master_task) + call broadcast_scalar(maxits_nonlin, master_task) + call broadcast_scalar(precond, master_task) + call broadcast_scalar(dim_fgmres, master_task) + call broadcast_scalar(dim_pgmres, master_task) + call broadcast_scalar(maxits_fgmres, master_task) + call broadcast_scalar(maxits_pgmres, master_task) + call broadcast_scalar(monitor_nonlin, master_task) + call broadcast_scalar(monitor_fgmres, master_task) + call broadcast_scalar(monitor_pgmres, master_task) + call broadcast_scalar(ortho_type, master_task) + call broadcast_scalar(reltol_nonlin, master_task) + call broadcast_scalar(reltol_fgmres, master_task) + call broadcast_scalar(reltol_pgmres, master_task) + call broadcast_scalar(algo_nonlin, master_task) + call broadcast_scalar(fpfunc_andacc, master_task) + call broadcast_scalar(dim_andacc, master_task) + call broadcast_scalar(reltol_andacc, master_task) + call broadcast_scalar(damping_andacc, master_task) + call broadcast_scalar(start_andacc, master_task) + call broadcast_scalar(use_mean_vrel, master_task) + call broadcast_scalar(conduct, master_task) + call broadcast_scalar(R_ice, master_task) + call broadcast_scalar(R_pnd, master_task) + call broadcast_scalar(R_snw, master_task) + call broadcast_scalar(dT_mlt, master_task) + call broadcast_scalar(rsnw_mlt, master_task) + call broadcast_scalar(kalg, master_task) + call broadcast_scalar(hp1, master_task) + call broadcast_scalar(hs0, master_task) + call broadcast_scalar(hs1, master_task) + call broadcast_scalar(dpscale, master_task) + call broadcast_scalar(frzpnd, master_task) + call broadcast_scalar(rfracmin, master_task) + call broadcast_scalar(rfracmax, master_task) + call broadcast_scalar(pndaspect, master_task) + call broadcast_scalar(albicev, master_task) + call broadcast_scalar(albicei, master_task) + call broadcast_scalar(albsnowv, master_task) + call broadcast_scalar(albsnowi, master_task) + call broadcast_scalar(ahmax, master_task) + call broadcast_scalar(atmbndy, master_task) + call broadcast_scalar(fyear_init, master_task) + call broadcast_scalar(ycycle, master_task) + call broadcast_scalar(atm_data_format, master_task) + call broadcast_scalar(atm_data_type, master_task) + call broadcast_scalar(atm_data_dir, master_task) + call broadcast_scalar(rotate_wind, master_task) + call broadcast_scalar(calc_strair, master_task) + call broadcast_scalar(calc_Tsfc, master_task) + call broadcast_scalar(formdrag, master_task) + call broadcast_scalar(highfreq, master_task) + call broadcast_scalar(natmiter, master_task) + call broadcast_scalar(atmiter_conv, master_task) + call broadcast_scalar(update_ocn_f, master_task) + call broadcast_scalar(l_mpond_fresh, master_task) + call broadcast_scalar(ustar_min, master_task) + call broadcast_scalar(emissivity, master_task) + call broadcast_scalar(fbot_xfer_type, master_task) + call broadcast_scalar(precip_units, master_task) + call broadcast_scalar(oceanmixed_ice, master_task) + call broadcast_scalar(wave_spec_type, master_task) + call broadcast_scalar(wave_spec_file, master_task) + call broadcast_scalar(nfreq, master_task) + call broadcast_scalar(tfrz_option, master_task) + call broadcast_scalar(ocn_data_format, master_task) + call broadcast_scalar(bgc_data_type, master_task) + call broadcast_scalar(fe_data_type, master_task) + call broadcast_scalar(ice_data_type, master_task) + call broadcast_scalar(bgc_data_dir, master_task) + call broadcast_scalar(ocn_data_type, master_task) + call broadcast_scalar(ocn_data_dir, master_task) + call broadcast_scalar(oceanmixed_file, master_task) + call broadcast_scalar(restore_ocn, master_task) + call broadcast_scalar(trestore, master_task) + call broadcast_scalar(restore_ice, master_task) + call broadcast_scalar(dbug, master_task) + call broadcast_array (latpnt(1:2), master_task) + call broadcast_array (lonpnt(1:2), master_task) + call broadcast_scalar(runid, master_task) + call broadcast_scalar(runtype, master_task) if (dbug) & ! else only master_task writes to file - call broadcast_scalar(nu_diag, master_task) + call broadcast_scalar(nu_diag, master_task) ! tracers - call broadcast_scalar(tr_iage, master_task) - call broadcast_scalar(restart_age, master_task) - call broadcast_scalar(tr_FY, master_task) - call broadcast_scalar(restart_FY, master_task) - call broadcast_scalar(tr_lvl, master_task) - call broadcast_scalar(restart_lvl, master_task) - call broadcast_scalar(tr_pond_cesm, master_task) - call broadcast_scalar(restart_pond_cesm, master_task) - call broadcast_scalar(tr_pond_lvl, master_task) - call broadcast_scalar(restart_pond_lvl, master_task) - call broadcast_scalar(tr_pond_topo, master_task) - call broadcast_scalar(restart_pond_topo, master_task) - call broadcast_scalar(tr_iso, master_task) - call broadcast_scalar(restart_iso, master_task) - call broadcast_scalar(tr_aero, master_task) - call broadcast_scalar(restart_aero, master_task) - call broadcast_scalar(tr_fsd, master_task) - call broadcast_scalar(restart_fsd, master_task) - call broadcast_scalar(ncat, master_task) - call broadcast_scalar(nfsd, master_task) - call broadcast_scalar(nilyr, master_task) - call broadcast_scalar(nslyr, master_task) - call broadcast_scalar(nblyr, master_task) - call broadcast_scalar(n_iso, master_task) - call broadcast_scalar(n_aero, master_task) - call broadcast_scalar(n_zaero, master_task) - call broadcast_scalar(n_algae, master_task) - call broadcast_scalar(n_doc, master_task) - call broadcast_scalar(n_dic, master_task) - call broadcast_scalar(n_don, master_task) - call broadcast_scalar(n_fed, master_task) - call broadcast_scalar(n_fep, master_task) - call broadcast_scalar(a_rapid_mode, master_task) - call broadcast_scalar(Rac_rapid_mode, master_task) - call broadcast_scalar(aspect_rapid_mode, master_task) - call broadcast_scalar(dSdt_slow_mode, master_task) - call broadcast_scalar(phi_c_slow_mode, master_task) - call broadcast_scalar(phi_i_mushy, master_task) - call broadcast_scalar(sw_redist, master_task) - call broadcast_scalar(sw_frac, master_task) - call broadcast_scalar(sw_dtemp, master_task) + call broadcast_scalar(tr_iage, master_task) + call broadcast_scalar(restart_age, master_task) + call broadcast_scalar(tr_FY, master_task) + call broadcast_scalar(restart_FY, master_task) + call broadcast_scalar(tr_lvl, master_task) + call broadcast_scalar(restart_lvl, master_task) + call broadcast_scalar(tr_pond_cesm, master_task) + call broadcast_scalar(restart_pond_cesm, master_task) + call broadcast_scalar(tr_pond_lvl, master_task) + call broadcast_scalar(restart_pond_lvl, master_task) + call broadcast_scalar(tr_pond_topo, master_task) + call broadcast_scalar(restart_pond_topo, master_task) + call broadcast_scalar(tr_iso, master_task) + call broadcast_scalar(restart_iso, master_task) + call broadcast_scalar(tr_aero, master_task) + call broadcast_scalar(restart_aero, master_task) + call broadcast_scalar(tr_fsd, master_task) + call broadcast_scalar(restart_fsd, master_task) + call broadcast_scalar(ncat, master_task) + call broadcast_scalar(nfsd, master_task) + call broadcast_scalar(nilyr, master_task) + call broadcast_scalar(nslyr, master_task) + call broadcast_scalar(nblyr, master_task) + call broadcast_scalar(n_iso, master_task) + call broadcast_scalar(n_aero, master_task) + call broadcast_scalar(n_zaero, master_task) + call broadcast_scalar(n_algae, master_task) + call broadcast_scalar(n_doc, master_task) + call broadcast_scalar(n_dic, master_task) + call broadcast_scalar(n_don, master_task) + call broadcast_scalar(n_fed, master_task) + call broadcast_scalar(n_fep, master_task) + call broadcast_scalar(a_rapid_mode, master_task) + call broadcast_scalar(floediam, master_task) + call broadcast_scalar(hfrazilmin, master_task) + call broadcast_scalar(Rac_rapid_mode, master_task) + call broadcast_scalar(aspect_rapid_mode, master_task) + call broadcast_scalar(dSdt_slow_mode, master_task) + call broadcast_scalar(phi_c_slow_mode, master_task) + call broadcast_scalar(phi_i_mushy, master_task) + call broadcast_scalar(sw_redist, master_task) + call broadcast_scalar(sw_frac, master_task) + call broadcast_scalar(sw_dtemp, master_task) #ifdef CESMCOUPLED pointer_file = trim(pointer_file) // trim(inst_suffix) @@ -846,11 +855,7 @@ subroutine input_data abort_list = trim(abort_list)//":1" endif -!deprecate upwind if (advection /= 'remap' .and. advection /= 'upwind' .and. advection /= 'none') then - if (advection /= 'remap' .and. advection /= 'none') then - if (trim(advection) == 'upwind') then - if (my_task == master_task) write(nu_diag,*) subname//' ERROR: upwind advection has been deprecated' - 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" endif @@ -889,6 +894,16 @@ subroutine input_data abort_list = trim(abort_list)//":33" endif + if (seabed_stress) then + if (seabed_stress_method /= 'LKD' .and. seabed_stress_method /= 'probabilistic') then + if (my_task == master_task) then + write(nu_diag,*) subname//' ERROR: invalid seabed stress method' + write(nu_diag,*) subname//' ERROR: seabed_stress_method should be LKD or probabilistic' + endif + abort_list = trim(abort_list)//":34" + endif + endif + rpcesm = 0 rplvl = 0 rptopo = 0 @@ -1139,203 +1154,237 @@ subroutine input_data if (my_task == master_task) then write(nu_diag,*) ' Overview of model configuration with relevant parameters' - write(nu_diag,*) ' ========================================================' - write(nu_diag,*) ' For details, compare namelist output below with the' - write(nu_diag,*) ' Case Settings section in the model documentation.' + write(nu_diag,*) '=========================================================' + write(nu_diag,*) 'For details, compare namelist output below with the' + write(nu_diag,*) 'Case Settings section in the model documentation.' write(nu_diag,*) ' ' write(nu_diag,*) ' Calendar' write(nu_diag,*) '--------------------------------' - write(nu_diag,1022) ' days_per_year = ',days_per_year,' number of days in a model year' + write(nu_diag,1020) ' days_per_year = ',days_per_year,' : number of days in a model year' if (use_leap_years) then - tmpstr2 = ' leap days are included' + tmpstr2 = ' : leap days are included' else - tmpstr2 = ' leap days are not included' + tmpstr2 = ' : leap days are not included' endif - write(nu_diag,1012) ' use_leap_years = ',use_leap_years,trim(tmpstr2) - write(nu_diag,1002) ' dt = ', dt, ' model time step' + write(nu_diag,1010) ' use_leap_years = ',use_leap_years,trim(tmpstr2) + write(nu_diag,1002) ' dt = ', dt, ' : model time step' write(nu_diag,*) ' ' write(nu_diag,*) ' Grid, Discretization' write(nu_diag,*) '--------------------------------' - if (trim(grid_type) == 'rectangular') & - write(nu_diag,*) 'grid_type = ', & - trim(grid_type),': internally defined, rectangular grid' - if (trim(grid_type) == 'regional') & - write(nu_diag,*) 'grid_type = ', & - trim(grid_type),': user-defined, regional grid' - if (trim(grid_type) == 'displaced_pole') & - write(nu_diag,*) 'grid_type = ', & - trim(grid_type),': user-defined grid with rotated north pole' - if (trim(grid_type) == 'tripole') then - write(nu_diag,*) 'grid_type = ', & - trim(grid_type),': user-defined grid with northern hemisphere zipper' - if (trim(ns_boundary_type) == 'tripole') then - tmpstr2 = ' on U points (nodes)' - elseif (trim(ns_boundary_type) == 'tripoleT') then - tmpstr2 = ' on T points (cell centers)' - endif - write(nu_diag,*) 'ns_boundary_type = ', trim(ns_boundary_type),trim(tmpstr2) - endif + tmpstr2 = ' ' + if (trim(grid_type) == 'rectangular') tmpstr2 = ' : internally defined, rectangular grid' + if (trim(grid_type) == 'regional') tmpstr2 = ' : user-defined, regional grid' + if (trim(grid_type) == 'displaced_pole') tmpstr2 = ' : user-defined grid with rotated north pole' + if (trim(grid_type) == 'tripole') tmpstr2 = ' : user-defined grid with northern hemisphere zipper' + write(nu_diag,1030) ' grid_type = ',trim(grid_type),trim(tmpstr2) if (trim(grid_type) /= 'rectangular') then if (use_bathymetry) then - tmpstr2 = ' bathymetric input data is used' + tmpstr2 = ' : bathymetric input data is used' else - tmpstr2 = ' bathymetric input data is not used' + tmpstr2 = ' : bathymetric input data is not used' endif - write(nu_diag,1012) ' use_bathymetry = ', use_bathymetry,trim(tmpstr2) - write(nu_diag,*) ' bathymetry_format= ', trim(bathymetry_format) + write(nu_diag,1010) ' use_bathymetry = ', use_bathymetry,trim(tmpstr2) + write(nu_diag,1030) ' bathymetry_format= ', trim(bathymetry_format) endif - write(nu_diag,1022) ' nilyr = ', nilyr, ' number of ice layers (equal thickness)' - write(nu_diag,1022) ' nslyr = ', nslyr, ' number of snow layers (equal thickness)' - write(nu_diag,1022) ' nblyr = ', nblyr, ' number of bio layers (equal thickness)' + write(nu_diag,1020) ' nilyr = ', nilyr, ' : number of ice layers (equal thickness)' + write(nu_diag,1020) ' nslyr = ', nslyr, ' : number of snow layers (equal thickness)' + write(nu_diag,1020) ' nblyr = ', nblyr, ' : number of bio layers (equal thickness)' if (trim(shortwave) == 'dEdd') & write(nu_diag,*) 'dEdd interior and sfc scattering layers are used in both ice, snow (unequal)' - write(nu_diag,1022) ' ncat = ', ncat, ' number of ice categories' + write(nu_diag,1020) ' ncat = ', ncat, ' : number of ice categories' if (kcatbound == 0) then - tmpstr2 = ' original ITD category bounds' + tmpstr2 = ' : original ITD category bounds' elseif (kcatbound == 1) then - tmpstr2 = ' round-number category bounds' + tmpstr2 = ' : round-number category bounds' elseif (kcatbound == 2) then - tmpstr2 = ' WMO standard ITD categories' + tmpstr2 = ' : WMO standard ITD categories' elseif (kcatbound == -1) then - tmpstr2 = ' one thickness category' + tmpstr2 = ' : one thickness category' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,1022) ' kcatbound = ', kcatbound,trim(tmpstr2) + write(nu_diag,1020) ' kcatbound = ', kcatbound,trim(tmpstr2) if (kitd==0) then - tmpstr2 = ' delta function ITD approx' + tmpstr2 = ' : delta function ITD approx' else - tmpstr2 = ' linear remapping ITD approx' + tmpstr2 = ' : linear remapping ITD approx' endif - write(nu_diag,1022) ' kitd = ', kitd,trim(tmpstr2) + write(nu_diag,1020) ' kitd = ', kitd,trim(tmpstr2) if (tr_fsd) then - tmpstr2 = ' floe size distribution is enabled' - ! write(nu_diag,1002) ' floediam = ', floediam, ' constant floe diameter' + tmpstr2 = ' : floe size distribution is enabled' else - tmpstr2 = ' floe size distribution is disabled' + tmpstr2 = ' : floe size distribution is disabled' endif - write(nu_diag,1012) ' tr_fsd = ', tr_fsd,trim(tmpstr2) - write(nu_diag,1022) ' nfsd = ', nfsd, ' number of floe size categories' + write(nu_diag,1010) ' tr_fsd = ', tr_fsd,trim(tmpstr2) + write(nu_diag,1020) ' nfsd = ', nfsd, ' : number of floe size categories' write(nu_diag,*) ' ' write(nu_diag,*) ' Horizontal Dynamics' write(nu_diag,*) '--------------------------------' if (kdyn == 1) then - tmpstr2 = ' elastic-viscous-plastic dynamics' + tmpstr2 = ' : elastic-viscous-plastic dynamics' elseif (kdyn == 2) then - tmpstr2 = ' elastic-anisotropic-plastic dynamics' + tmpstr2 = ' : elastic-anisotropic-plastic dynamics' elseif (kdyn == 3) then - tmpstr2 = ' viscous-plastic dynamics' + tmpstr2 = ' : viscous-plastic dynamics' elseif (kdyn < 1) then - tmpstr2 = ' dynamics disabled' + tmpstr2 = ' : dynamics disabled' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,1022) ' kdyn = ', kdyn,trim(tmpstr2) + write(nu_diag,1020) ' kdyn = ', kdyn,trim(tmpstr2) if (kdyn >= 1) then if (kdyn == 1 .or. kdyn == 2) then if (revised_evp) then - tmpstr2 = ' revised EVP formulation used' - write(nu_diag,1007) ' arlx = ', arlx, ' stress equation factor alpha' - write(nu_diag,1007) ' brlx = ', brlx, ' stress equation factor beta' + tmpstr2 = ' : revised EVP formulation used' + write(nu_diag,1002) ' arlx = ', arlx, ' : stress equation factor alpha' + write(nu_diag,1002) ' brlx = ', brlx, ' : stress equation factor beta' else - tmpstr2 = ' revised EVP formulation not used' + tmpstr2 = ' : revised EVP formulation not used' endif - write(nu_diag,1012) ' revised_evp = ', revised_evp,trim(tmpstr2) - write(nu_diag,1022) ' kevp_kernel = ', kevp_kernel,' EVP solver' + write(nu_diag,1010) ' revised_evp = ', revised_evp,trim(tmpstr2) - write(nu_diag,1022) ' ndtd = ', ndtd, ' number of dynamics/advection/ridging/steps per thermo timestep' - write(nu_diag,1022) ' ndte = ', ndte, ' number of EVP or EAP subcycles' - endif + if (kevp_kernel == 0) then + tmpstr2 = ' : original EVP solver' + elseif (kevp_kernel == 2 .or. kevp_kernel == 102) then + tmpstr2 = ' : vectorized EVP solver' + else + tmpstr2 = ' : unknown value' + endif + write(nu_diag,1020) ' kevp_kernel = ', kevp_kernel,trim(tmpstr2) + + write(nu_diag,1020) ' ndtd = ', ndtd, ' : number of dynamics/advection/ridging/steps per thermo timestep' + write(nu_diag,1020) ' ndte = ', ndte, ' : number of EVP or EAP subcycles' + endif - if (kdyn == 1 .or. kdyn == 3) then - write(nu_diag,*) 'yield_curve = ', trim(yield_curve) - if (trim(yield_curve) == 'ellipse') & - write(nu_diag,1007) ' e_ratio = ', e_ratio, ' aspect ratio of ellipse' - endif + 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' + endif if (trim(coriolis) == 'latitude') then - tmpstr2 = ': latitude-dependent Coriolis parameter' + tmpstr2 = ' : latitude-dependent Coriolis parameter' elseif (trim(coriolis) == 'contant') then - tmpstr2 = ' = 1.46e-4/s' + tmpstr2 = ' : = 1.46e-4/s' elseif (trim(coriolis) == 'zero') then - tmpstr2 = ' = 0.0' + tmpstr2 = ' : = 0.0' + else + tmpstr2 = ': unknown value' endif - write(nu_diag,*) 'coriolis = ',trim(coriolis),trim(tmpstr2) + write(nu_diag,1030) ' coriolis = ',trim(coriolis),trim(tmpstr2) if (trim(ssh_stress) == 'geostrophic') then - tmpstr2 = ': from ocean velocity' + tmpstr2 = ' : from ocean velocity' elseif (trim(ssh_stress) == 'coupled') then - tmpstr2 = ': from coupled sea surface height gradients' + tmpstr2 = ' : from coupled sea surface height gradients' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) 'ssh_stress = ',trim(ssh_stress),trim(tmpstr2) - - if (ktransport == 1) then - tmpstr2 = ' transport enabled' - if (trim(advection) == 'remap') then - tmpstr2 = ': linear remapping advection' -!deprecate upwind elseif (trim(advection) == 'upwind') then -!deprecate upwind tmpstr2 = ': donor cell (upwind) advection' - elseif (trim(advection) == 'none') then - tmpstr2 = ': advection off' - endif - write(nu_diag,*) 'advection = ', trim(advection),trim(tmpstr2) + 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 + tmpstr2 = ' : donor cell (upwind) advection' + elseif (trim(advection) == 'none') then + tmpstr2 = ' : advection disabled by ktransport namelist' else - tmpstr2 = ' transport disabled' + tmpstr2 = ' : unknown value' endif - write(nu_diag,1022) ' ktransport = ', ktransport,trim(tmpstr2) + write(nu_diag,1030) ' advection = ', trim(advection),trim(tmpstr2) - if (basalstress) then - tmpstr2 = ' use basal stress parameterization for landfast ice' + if (seabed_stress) then + tmpstr2 = ' : use seabed stress parameterization for landfast ice' else - tmpstr2 = ' basal stress not used for landfast ice' + tmpstr2 = ' : no seabed stress parameterization' endif - write(nu_diag,1012) ' basalstress = ', basalstress,trim(tmpstr2) - if (basalstress) then - write(nu_diag,1007) ' k1 = ', k1, ' free parameter for landfast ice' - write(nu_diag,1007) ' k2 = ', k2, ' free parameter for landfast ice' - write(nu_diag,1007) ' alphab = ', alphab, ' factor for landfast ice' - write(nu_diag,1007) ' threshold_hw = ', threshold_hw, ' max water depth for grounding ice' + write(nu_diag,1010) ' seabed_stress = ', seabed_stress,trim(tmpstr2) + if (seabed_stress) then + write(nu_diag,1030) ' seabed method = ',trim(seabed_stress_method) + if (seabed_stress_method == 'LKD') then + write(nu_diag,1002) ' k1 = ', k1, ' : free parameter for landfast ice' + write(nu_diag,1002) ' k2 = ', k2, ' : free parameter for landfast ice' + write(nu_diag,1002) ' alphab = ', alphab, ' : factor for landfast ice' + write(nu_diag,1002) ' threshold_hw = ', threshold_hw, ' : max water depth for grounding ice' + elseif (seabed_stress_method == 'probabilistic') then + write(nu_diag,1002) ' alphab = ', alphab, ' : factor for landfast ice' + endif endif - write(nu_diag,1007) ' Ktens = ', Ktens, ' tensile strength factor' + write(nu_diag,1002) ' Ktens = ', Ktens, ' : tensile strength factor' + + if (kdyn == 3) then + write(nu_diag,1020) ' maxits_nonlin = ', maxits_nonlin,' : max nb of iteration for nonlinear solver' + write(nu_diag,1030) ' precond = ', trim(precond),' : preconditioner for FGMRES' + write(nu_diag,1020) ' dim_fgmres = ', dim_fgmres,' : size of FGMRES Krylov subspace' + write(nu_diag,1020) ' dim_pgmres = ', dim_pgmres,' : size of PGMRES Krylov subspace' + write(nu_diag,1020) ' maxits_fgmres = ', maxits_fgmres,' : max nb of iteration for FGMRES' + write(nu_diag,1020) ' maxits_pgmres = ', maxits_pgmres,' : max nb of iteration for PGMRES' + write(nu_diag,1010) ' monitor_nonlin = ', monitor_nonlin,' : print nonlinear residual norm' + write(nu_diag,1010) ' monitor_fgmres = ', monitor_fgmres,' : print FGMRES residual norm' + write(nu_diag,1010) ' monitor_pgmres = ', monitor_pgmres,' : print PGMRES residual norm' + write(nu_diag,1030) ' ortho_type = ', trim(ortho_type),' : type of orthogonalization for FGMRES' + write(nu_diag,1009) ' reltol_nonlin = ', reltol_nonlin,' : nonlinear stopping criterion' + write(nu_diag,1009) ' reltol_fgmres = ', reltol_fgmres,' : FGMRES stopping criterion' + write(nu_diag,1009) ' reltol_pgmres = ', reltol_pgmres,' : PGMRES stopping criterion' + write(nu_diag,1030) ' algo_nonlin = ', trim(algo_nonlin),' : nonlinear algorithm' + write(nu_diag,1010) ' use_mean_vrel = ', use_mean_vrel,' : use mean of previous 2 iterates to compute vrel' + if (algo_nonlin == 'anderson') then + write(nu_diag,1020) ' fpfunc_andacc = ', fpfunc_andacc,' : fixed point function for Anderson acceleration' + write(nu_diag,1020) ' dim_andacc = ', dim_andacc,' : size of Anderson minimization matrix' + write(nu_diag,1009) ' reltol_andacc = ', reltol_andacc,' : relative tolerance for Anderson acceleration' + write(nu_diag,1000) ' damping_andacc = ', damping_andacc,' : damping factor for Anderson acceleration' + write(nu_diag,1020) ' start_andacc = ', start_andacc,' : nonlinear iteration at which acceleration starts' + endif + endif + endif ! kdyn enabled write(nu_diag,*) ' ' write(nu_diag,*) ' Mechanical Deformation (Ridging) and Ice Strength' write(nu_diag,*) '--------------------------------------------------' if (kridge == 1) then - tmpstr2 = ' ridging enabled' + tmpstr2 = ' : ridging enabled' else - tmpstr2 = ' ridging disabled' + tmpstr2 = ' : ridging disabled' endif - write(nu_diag,1012) ' tr_lvl = ', tr_lvl,' ridging related tracers' - write(nu_diag,1022) ' kridge = ', kridge,trim(tmpstr2) + write(nu_diag,1010) ' tr_lvl = ', tr_lvl,' : ridging related tracers' + write(nu_diag,1020) ' kridge = ', kridge,trim(tmpstr2) if (kridge == 1) then if (krdg_partic == 1) then - tmpstr2 = ' new participation function' + tmpstr2 = ' : new participation function' else - tmpstr2 = ' old participation function' + tmpstr2 = ' : old participation function' endif - write(nu_diag,1022) ' krdg_partic = ', krdg_partic,trim(tmpstr2) + write(nu_diag,1020) ' krdg_partic = ', krdg_partic,trim(tmpstr2) if (krdg_partic == 1) & - write(nu_diag,1007) ' mu_rdg = ', mu_rdg,' e-folding scale of ridged ice' + write(nu_diag,1002) ' mu_rdg = ', mu_rdg,' : e-folding scale of ridged ice' if (krdg_redist == 1) then - tmpstr2 = ' new redistribution function' + tmpstr2 = ' : new redistribution function' else - tmpstr2 = ' old redistribution function' + tmpstr2 = ' : old redistribution function' endif - write(nu_diag,1022) ' krdg_redist = ', krdg_redist,trim(tmpstr2) + write(nu_diag,1020) ' krdg_redist = ', krdg_redist,trim(tmpstr2) endif if (kstrength == 0) then - tmpstr2 = ' Hibler (1979)' + tmpstr2 = ' : Hibler (1979)' elseif (kstrength == 1) then - tmpstr2 = ' Rothrock (1975)' + tmpstr2 = ' : Rothrock (1975)' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,1022) ' kstrength = ', kstrength,trim(tmpstr2) + write(nu_diag,1020) ' kstrength = ', kstrength,trim(tmpstr2) if (kstrength == 0) then - write(nu_diag,1009) ' Pstar = ', Pstar, ' P* strength factor' - write(nu_diag,1007) ' Cstar = ', Cstar, ' C* strength exponent factor' + write(nu_diag,1009) ' Pstar = ', Pstar, ' : P* strength factor' + write(nu_diag,1002) ' Cstar = ', Cstar, ' : C* strength exponent factor' elseif (kstrength == 1) then - write(nu_diag,1007) ' Cf = ', Cf, ' ratio of ridging work to PE change' + write(nu_diag,1002) ' Cf = ', Cf, ' : ratio of ridging work to PE change' endif write(nu_diag,*) ' ' @@ -1343,186 +1392,206 @@ subroutine input_data write(nu_diag,*) '--------------------------------' if (ktherm == 1) then - tmpstr2 = ' Bitz and Lipscomb 1999 thermo' + tmpstr2 = ' : Bitz and Lipscomb 1999 thermo' elseif (ktherm == 2) then - tmpstr2 = ' mushy-layer thermo' + tmpstr2 = ' : mushy-layer thermo' elseif (ktherm == 0) then - tmpstr2 = ' zero-layer thermo' + tmpstr2 = ' : zero-layer thermo' elseif (ktherm < 0) then - tmpstr2 = ' thermodynamics disabled' + tmpstr2 = ' : Thermodynamics disabled' + else + tmpstr2 = ' : unknown value' endif + write(nu_diag,1020) ' ktherm = ', ktherm,trim(tmpstr2) if (ktherm >= 0) then - write(nu_diag,1022) ' ktherm = ', ktherm,trim(tmpstr2) - write(nu_diag,1002) ' dt = ', dt, ' thermodynamic time step' - write(nu_diag,1007) ' ksno = ', ksno,' snow thermal conductivity' + write(nu_diag,1002) ' dt = ', dt, ' : thermodynamic time step' + write(nu_diag,1002) ' ksno = ', ksno,' : snow thermal conductivity' if (ktherm == 1) & - write(nu_diag,*) 'conduct = ', trim(conduct),' ice thermal conductivity' - write(nu_diag,1012) ' sw_redist = ', sw_redist,' redistribute internal shortwave to surface' - write(nu_diag,1002) ' sw_frac = ', sw_frac,' fraction redistributed' - write(nu_diag,1002) ' sw_dtemp = ', sw_dtemp,' temperature difference from freezing to redistribute' + write(nu_diag,1030) ' conduct = ', trim(conduct),' : ice thermal conductivity' if (ktherm == 2) then - write(nu_diag,1002) ' a_rapid_mode = ', a_rapid_mode,' brine channel diameter' - write(nu_diag,1007) ' Rac_rapid_mode = ', Rac_rapid_mode,' critical Rayleigh number' - write(nu_diag,1007) ' aspect_rapid_mode= ', aspect_rapid_mode,' brine convection aspect ratio' - write(nu_diag,*) 'dSdt_slow_mode = ', dSdt_slow_mode,' drainage strength parameter' - write(nu_diag,1007) ' phi_c_slow_mode = ', phi_c_slow_mode,' critical liquid fraction' - write(nu_diag,1007) ' phi_i_mushy = ', phi_i_mushy,' solid fraction at lower boundary' + write(nu_diag,1002) ' a_rapid_mode = ', a_rapid_mode,' : brine channel diameter' + write(nu_diag,1002) ' Rac_rapid_mode = ', Rac_rapid_mode,' : critical Rayleigh number' + write(nu_diag,1002) ' aspect_rapid_mode= ', aspect_rapid_mode,' : brine convection aspect ratio' + write(nu_diag,1009) ' dSdt_slow_mode = ', dSdt_slow_mode,' : drainage strength parameter' + write(nu_diag,1002) ' phi_c_slow_mode = ', phi_c_slow_mode,' : critical liquid fraction' + write(nu_diag,1002) ' phi_i_mushy = ', phi_i_mushy,' : solid fraction at lower boundary' endif - endif - !write(nu_diag,1007) ' hfrazilmin = ', hfrazilmin,' minimum new frazil ice thickness' - - write(nu_diag,*) ' ' - write(nu_diag,*) ' Radiation' - write(nu_diag,*) '--------------------------------' - if (trim(shortwave) == 'dEdd') then - tmpstr2 = ': delta-Eddington multiple-scattering method' - elseif (trim(shortwave) == 'ccsm3') then - tmpstr2 = ': NCAR CCSM3 distribution method' - endif - write(nu_diag,*) ' shortwave = ', trim(shortwave),trim(tmpstr2) - if (trim(shortwave) == 'dEdd') then - write(nu_diag,1007) ' R_ice = ', R_ice,' tuning parameter for sea ice albedo' - write(nu_diag,1007) ' R_pnd = ', R_pnd,' tuning parameter for ponded sea ice albedo' - write(nu_diag,1007) ' R_snw = ', R_snw,' tuning parameter for snow broadband albedo' - write(nu_diag,1007) ' dT_mlt = ', dT_mlt,' change in temperature per change in snow grain radius' - write(nu_diag,1002) ' rsnw_mlt = ', rsnw_mlt,' maximum melting snow grain radius' - write(nu_diag,1007) ' kalg = ', kalg,' absorption coefficient for algae' - else - if (trim(albedo_type) == 'ccsm3') then - tmpstr2 = ': NCAR CCSM3 albedos' - elseif (trim(albedo_type) == 'constant') then - tmpstr2 = ': four constant albedos' + write(nu_diag,1002) ' hfrazilmin = ', hfrazilmin,' : minimum new frazil ice thickness' + + write(nu_diag,*) ' ' + write(nu_diag,*) ' Radiation' + write(nu_diag,*) '--------------------------------' + if (trim(shortwave) == 'dEdd') then + tmpstr2 = ' : delta-Eddington multiple-scattering method' + elseif (trim(shortwave) == 'ccsm3') then + tmpstr2 = ' : NCAR CCSM3 distribution method' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) 'albedo_type = ', trim(albedo_type),trim(tmpstr2) - if (trim(albedo_type) == 'ccsm3') then - write(nu_diag,1007) ' albicev = ', albicev,' visible ice albedo for thicker ice' - write(nu_diag,1007) ' albicei = ', albicei,' near infrared ice albedo for thicker ice' - write(nu_diag,1007) ' albsnowv = ', albsnowv,' visible, cold snow albedo' - write(nu_diag,1007) ' albsnowi = ', albsnowi,' near infrared, cold snow albedo' - write(nu_diag,1007) ' ahmax = ', ahmax,' albedo is constant above this thickness' + write(nu_diag,1030) ' shortwave = ', trim(shortwave),trim(tmpstr2) + if (trim(shortwave) == 'dEdd') then + write(nu_diag,1002) ' R_ice = ', R_ice,' : tuning parameter for sea ice albedo' + write(nu_diag,1002) ' R_pnd = ', R_pnd,' : tuning parameter for ponded sea ice albedo' + write(nu_diag,1002) ' R_snw = ', R_snw,' : tuning parameter for snow broadband albedo' + write(nu_diag,1002) ' dT_mlt = ', dT_mlt,' : change in temperature per change in snow grain radius' + write(nu_diag,1002) ' rsnw_mlt = ', rsnw_mlt,' : maximum melting snow grain radius' + write(nu_diag,1002) ' kalg = ', kalg,' : absorption coefficient for algae' + else + if (trim(albedo_type) == 'ccsm3') then + tmpstr2 = ' : NCAR CCSM3 albedos' + elseif (trim(albedo_type) == 'constant') then + tmpstr2 = ' : four constant albedos' + else + tmpstr2 = ' : unknown value' + endif + write(nu_diag,1030) ' albedo_type = ', trim(albedo_type),trim(tmpstr2) + if (trim(albedo_type) == 'ccsm3') then + write(nu_diag,1002) ' albicev = ', albicev,' : visible ice albedo for thicker ice' + write(nu_diag,1002) ' albicei = ', albicei,' : near infrared ice albedo for thicker ice' + write(nu_diag,1002) ' albsnowv = ', albsnowv,' : visible, cold snow albedo' + write(nu_diag,1002) ' albsnowi = ', albsnowi,' : near infrared, cold snow albedo' + write(nu_diag,1002) ' ahmax = ', ahmax,' : albedo is constant above this thickness' + write(nu_diag,1002) ' ahmax = ', ahmax,' : albedo is constant above this thickness' + endif + endif + write(nu_diag,1000) ' emissivity = ', emissivity,' : emissivity of snow and ice' + write(nu_diag,1010) ' sw_redist = ', sw_redist,' : redistribute internal shortwave to surface' + if (sw_redist) then + write(nu_diag,1002) ' sw_frac = ', sw_frac,' : fraction redistributed' + write(nu_diag,1002) ' sw_dtemp = ', sw_dtemp,' : temperature difference from freezing to redistribute' endif endif - write(nu_diag,1007) ' emissivity = ', emissivity,' emissivity of snow and ice' write(nu_diag,*) ' ' write(nu_diag,*) ' Atmospheric Forcing / Coupling' write(nu_diag,*) '--------------------------------' - write(nu_diag,1012) ' calc_Tsfc = ', calc_Tsfc,' calculate surface temperature as part of thermo' - write(nu_diag,1012) ' calc_strair = ', calc_strair,' calculate wind stress and speed' - write(nu_diag,1012) ' rotate_wind = ', rotate_wind,' rotate wind/stress to computational grid' - write(nu_diag,1012) ' formdrag = ', formdrag,' use form drag parameterization' + write(nu_diag,1010) ' calc_Tsfc = ', calc_Tsfc,' : calculate surface temperature as part of thermo' + write(nu_diag,1010) ' calc_strair = ', calc_strair,' : calculate wind stress and speed' + write(nu_diag,1010) ' rotate_wind = ', rotate_wind,' : rotate wind/stress to computational grid' + write(nu_diag,1010) ' formdrag = ', formdrag,' : use form drag parameterization' if (trim(atmbndy) == 'default') then - tmpstr2 = ': stability-based boundary layer' - write(nu_diag,1012) ' highfreq = ', highfreq,' high-frequency atmospheric coupling' - write(nu_diag,1022) ' natmiter = ', natmiter,' number of atmo boundary layer iterations' - write(nu_diag,1006) ' atmiter_conv = ', atmiter_conv,' convergence criterion for ustar' + tmpstr2 = ' : stability-based boundary layer' + 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' + tmpstr2 = ' : boundary layer uses bulk transfer coefficients' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) 'atmbndy = ', trim(atmbndy),trim(tmpstr2) + write(nu_diag,1030) ' atmbndy = ', trim(atmbndy),trim(tmpstr2) write(nu_diag,*) ' ' write(nu_diag,*) ' Oceanic Forcing / Coupling' write(nu_diag,*) '--------------------------------' if (oceanmixed_ice) then - tmpstr2 = ' ocean mixed layer calculation (SST) enabled' + tmpstr2 = ' : ocean mixed layer calculation (SST) enabled' else - tmpstr2 = ' ocean mixed layer calculation (SST) disabled' + tmpstr2 = ' : ocean mixed layer calculation (SST) disabled' endif - write(nu_diag,1012) ' oceanmixed_ice = ', oceanmixed_ice,trim(tmpstr2) + write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice,trim(tmpstr2) if (oceanmixed_ice) then write(nu_diag,*) ' WARNING: ocean mixed layer ON' write(nu_diag,*) ' WARNING: will impact ocean forcing interaction' write(nu_diag,*) ' WARNING: coupled forcing will be modified by mixed layer routine' endif if (trim(tfrz_option) == 'minus1p8') then - tmpstr2 = ': constant ocean freezing temperature (-1.8C)' + tmpstr2 = ' : constant ocean freezing temperature (-1.8C)' elseif (trim(tfrz_option) == 'linear_salt') then - tmpstr2 = ': linear function of salinity (use with ktherm=1)' + tmpstr2 = ' : linear function of salinity (use with ktherm=1)' elseif (trim(tfrz_option) == 'mushy') then - tmpstr2 = ': Assur (1958) as in mushy-layer thermo (ktherm=2)' + tmpstr2 = ' : Assur (1958) as in mushy-layer thermo (ktherm=2)' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) 'tfrz_option = ', trim(tfrz_option),trim(tmpstr2) + write(nu_diag,1030) ' tfrz_option = ', trim(tfrz_option),trim(tmpstr2) if (update_ocn_f) then - tmpstr2 = ' frazil water/salt fluxes included in ocean fluxes' + tmpstr2 = ' : frazil water/salt fluxes included in ocean fluxes' else - tmpstr2 = ' frazil water/salt fluxes not included in ocean fluxes' + tmpstr2 = ' : frazil water/salt fluxes not included in ocean fluxes' endif - write(nu_diag,1012) ' update_ocn_f = ', update_ocn_f,trim(tmpstr2) + write(nu_diag,1010) ' update_ocn_f = ', update_ocn_f,trim(tmpstr2) if (l_mpond_fresh .and. tr_pond_topo) then - tmpstr2 = ' retain (topo) pond water until ponds drain' + tmpstr2 = ' : retain (topo) pond water until ponds drain' else - tmpstr2 = ' pond water not retained on ice (virtual only)' + tmpstr2 = ' : pond water not retained on ice (virtual only)' endif - write(nu_diag,1012) ' l_mpond_fresh = ', l_mpond_fresh,trim(tmpstr2) + write(nu_diag,1010) ' l_mpond_fresh = ', l_mpond_fresh,trim(tmpstr2) if (trim(fbot_xfer_type) == 'constant') then - tmpstr2 = ': ocean heat transfer coefficient is constant' + tmpstr2 = ' : ocean heat transfer coefficient is constant' elseif (trim(fbot_xfer_type) == 'Cdn_ocn') then - tmpstr2 = ': variable ocean heat transfer coefficient' ! only used with form_drag=T? + tmpstr2 = ' : variable ocean heat transfer coefficient' ! only used with form_drag=T? + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) 'fbot_xfer_type = ', trim(fbot_xfer_type),trim(tmpstr2) - write(nu_diag,1006) ' ustar_min = ', ustar_min,' minimum value of ocean friction velocity' + write(nu_diag,1030) ' fbot_xfer_type = ', trim(fbot_xfer_type),trim(tmpstr2) + write(nu_diag,1000) ' ustar_min = ', ustar_min,' : minimum value of ocean friction velocity' if (tr_fsd) then - if (wave_spec) then - tmpstr2 = ' use wave spectrum for floe size distribution' - else - tmpstr2 = ' floe size distribution does not use wave spectrum' - endif - write(nu_diag,1012) ' wave_spec = ', wave_spec,trim(tmpstr2) - if (wave_spec) then - if (trim(wave_spec_type) == 'none') then - tmpstr2 = ': no wave data provided, no wave-ice interactions' - elseif (trim(wave_spec_type) == 'profile') then - tmpstr2 = ': use fixed dummy wave spectrum for testing' - elseif (trim(wave_spec_type) == 'constant') then - tmpstr2 = ': constant wave spectrum data file provided for testing' - elseif (trim(wave_spec_type) == 'random') then - tmpstr2 = ': wave data file provided, spectrum generated using random number' + write(nu_diag,1002) ' floediam = ', floediam, ' constant floe diameter' + if (wave_spec) then + tmpstr2 = ' : use wave spectrum for floe size distribution' + else + tmpstr2 = ' : floe size distribution does not use wave spectrum' endif - write(nu_diag,*) 'wave_spec_type = ', trim(wave_spec_type),trim(tmpstr2) - endif - write(nu_diag,1022) ' nfreq = ', nfreq,' number of wave spectral forcing frequencies' + write(nu_diag,1010) ' wave_spec = ', wave_spec,trim(tmpstr2) + if (wave_spec) then + if (trim(wave_spec_type) == 'none') then + tmpstr2 = ' : no wave data provided, no wave-ice interactions' + elseif (trim(wave_spec_type) == 'profile') then + tmpstr2 = ' : use fixed dummy wave spectrum for testing' + elseif (trim(wave_spec_type) == 'constant') then + tmpstr2 = ' : constant wave spectrum data file provided for testing' + elseif (trim(wave_spec_type) == 'random') then + tmpstr2 = ' : wave data file provided, spectrum generated using random number' + else + tmpstr2 = ' : unknown value' + endif + write(nu_diag,1030) ' wave_spec_type = ', trim(wave_spec_type),trim(tmpstr2) + endif + write(nu_diag,1020) ' nfreq = ', nfreq,' : number of wave spectral forcing frequencies' endif write(nu_diag,*) ' ' write(nu_diag,*) ' Age related tracers' write(nu_diag,*) '--------------------------------' - write(nu_diag,1012) ' tr_iage = ', tr_iage,' chronological ice age' - write(nu_diag,1012) ' tr_FY = ', tr_FY,' first-year ice area' + write(nu_diag,1010) ' tr_iage = ', tr_iage,' : chronological ice age' + write(nu_diag,1010) ' tr_FY = ', tr_FY,' : first-year ice area' write(nu_diag,*) ' ' write(nu_diag,*) ' Melt ponds' write(nu_diag,*) '--------------------------------' if (tr_pond_cesm) then - write(nu_diag,1012) ' tr_pond_cesm = ', tr_pond_cesm,' CESM pond formulation' - write(nu_diag,1007) ' pndaspect = ', pndaspect + write(nu_diag,1010) ' tr_pond_cesm = ', tr_pond_cesm,' : CESM pond formulation' + write(nu_diag,1002) ' pndaspect = ', pndaspect,' : ratio of pond depth to area fraction' elseif (tr_pond_lvl) then - write(nu_diag,1012) ' tr_pond_lvl = ', tr_pond_lvl,' level-ice pond formulation' - write(nu_diag,1007) ' pndaspect = ', pndaspect - write(nu_diag,1006) ' dpscale = ', dpscale,' time scale for flushing in permeable ice' + write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl,' : level-ice pond formulation' + write(nu_diag,1002) ' pndaspect = ', pndaspect,' : ratio of pond depth to area fraction' + write(nu_diag,1000) ' dpscale = ', dpscale,' : time scale for flushing in permeable ice' if (trim(frzpnd) == 'hlid') then - tmpstr2 = ': Stefan refreezing with pond ice thickness' + tmpstr2 = ' : Stefan refreezing with pond ice thickness' elseif (trim(frzpnd) == 'cesm') then - tmpstr2 = ': CESM refreezing empirical formula' + tmpstr2 = ' : CESM refreezing empirical formula' + else + tmpstr2 = ' : unknown value' endif - write(nu_diag,*) ' frzpnd = ', trim(frzpnd),trim(tmpstr2) - write(nu_diag,1007) ' hs1 = ', hs1,' snow depth of transition to pond ice' + write(nu_diag,1030) ' frzpnd = ', trim(frzpnd),trim(tmpstr2) + write(nu_diag,1002) ' hs1 = ', hs1,' : snow depth of transition to pond ice' elseif (tr_pond_topo) then - write(nu_diag,1012) ' tr_pond_topo = ', tr_pond_topo,' topo pond formulation' - write(nu_diag,1007) ' hp1 = ', hp1,' critical ice lid thickness for topo ponds' + write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo,' : topo pond formulation' + write(nu_diag,1002) ' hp1 = ', hp1,' : critical ice lid thickness for topo ponds' elseif (trim(shortwave) == 'ccsm3') then write(nu_diag,*) 'Pond effects on radiation are treated implicitly in the ccsm3 shortwave scheme' else - write(nu_diag,*) ' Using default dEdd melt pond scheme for testing only' + write(nu_diag,*) 'Using default dEdd melt pond scheme for testing only' endif if (trim(shortwave) == 'dEdd') then - write(nu_diag,1007) ' hs0 = ', hs0,' snow depth of transition to bare sea ice' + write(nu_diag,1002) ' hs0 = ', hs0,' : snow depth of transition to bare sea ice' endif - write(nu_diag,1007) ' rfracmin = ', rfracmin,' minimum fraction of melt water added to ponds' - write(nu_diag,1007) ' rfracmax = ', rfracmax,' maximum fraction of melt water added to ponds' + write(nu_diag,1002) ' rfracmin = ', rfracmin,' : minimum fraction of melt water added to ponds' + write(nu_diag,1002) ' rfracmax = ', rfracmax,' : maximum fraction of melt water added to ponds' write(nu_diag,*) ' ' write(nu_diag,*) ' Primary state variables, tracers' @@ -1531,191 +1600,135 @@ subroutine input_data write(nu_diag,*) 'Conserved properties (all tracers are conserved):' write(nu_diag,*) 'ice concentration, volume and enthalpy' write(nu_diag,*) 'snow volume and enthalpy' - if (ktherm == 2) write(nu_diag,*) 'ice salinity' - if (tr_fsd) write(nu_diag,1012) ' tr_fsd = ', tr_fsd,' floe size distribution' - if (tr_lvl) write(nu_diag,1012) ' tr_lvl = ', tr_lvl,' ridging related tracers' - if (tr_pond_lvl) write(nu_diag,1012) ' tr_pond_lvl = ', tr_pond_lvl,' level-ice pond formulation' - if (tr_pond_topo) write(nu_diag,1012) ' tr_pond_topo = ', tr_pond_topo,' topo pond formulation' - if (tr_pond_cesm) write(nu_diag,1012) ' tr_pond_cesm = ', tr_pond_cesm,' CESM pond formulation' - if (tr_iage) write(nu_diag,1012) ' tr_iage = ', tr_iage,' chronological ice age' - if (tr_FY) write(nu_diag,1012) ' tr_FY = ', tr_FY,' first-year ice area' - if (tr_iso) write(nu_diag,1012) ' tr_iso = ', tr_iso,' diagnostic isotope tracers' - if (tr_aero) write(nu_diag,1012) ' tr_aero = ', tr_aero,' CESM aerosol tracers' + if (ktherm == 2) write(nu_diag,1030) ' ice salinity' + if (tr_fsd) write(nu_diag,1010) ' tr_fsd = ', tr_fsd,' : floe size distribution' + if (tr_lvl) write(nu_diag,1010) ' tr_lvl = ', tr_lvl,' : ridging related tracers' + if (tr_pond_lvl) write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl,' : level-ice pond formulation' + if (tr_pond_topo) write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo,' : topo pond formulation' + if (tr_pond_cesm) write(nu_diag,1010) ' tr_pond_cesm = ', tr_pond_cesm,' : CESM pond formulation' + if (tr_iage) write(nu_diag,1010) ' tr_iage = ', tr_iage,' : chronological ice age' + if (tr_FY) write(nu_diag,1010) ' tr_FY = ', tr_FY,' : first-year ice area' + if (tr_iso) write(nu_diag,1010) ' tr_iso = ', tr_iso,' : diagnostic isotope tracers' + if (tr_aero) write(nu_diag,1010) ' tr_aero = ', tr_aero,' : CESM aerosol tracers' write(nu_diag,*) 'Non-conserved properties:' write(nu_diag,*) 'ice surface temperature' write(nu_diag,*) 'ice velocity components and internal stress' write(nu_diag,*) ' ' write(nu_diag,*) ' Other ice_in namelist parameters:' - write(nu_diag,*) ' ==================================== ' - write(nu_diag,*) ' ' + write(nu_diag,*) '===================================== ' if (trim(runid) /= 'unknown') & - write(nu_diag,*) ' runid = ', & - trim(runid) - write(nu_diag,1030) ' runtype = ', & - trim(runtype) - write(nu_diag,1020) ' year_init = ', year_init - write(nu_diag,1020) ' istep0 = ', istep0 - write(nu_diag,1020) ' npt = ', npt - write(nu_diag,1020) ' diagfreq = ', diagfreq - write(nu_diag,1010) ' print_global = ', print_global - write(nu_diag,1010) ' print_points = ', print_points - write(nu_diag,1030) ' bfbflag = ', bfbflag - write(nu_diag,1020) ' numin = ', numin - write(nu_diag,1020) ' numax = ', numax - write(nu_diag,1050) ' histfreq = ', histfreq(:) - write(nu_diag,1040) ' histfreq_n = ', histfreq_n(:) - write(nu_diag,1010) ' hist_avg = ', hist_avg - if (.not. hist_avg) write(nu_diag,*) ' History data will be snapshots' - write(nu_diag,*) ' history_dir = ', & - trim(history_dir) - write(nu_diag,*) ' history_file = ', & - trim(history_file) - write(nu_diag,1020) ' history_precision = ', history_precision - write(nu_diag,*) ' history_format = ', & - trim(history_format) + write(nu_diag,1031) ' runid = ', trim(runid) + write(nu_diag,1031) ' runtype = ', trim(runtype) + write(nu_diag,1021) ' year_init = ', year_init + write(nu_diag,1021) ' istep0 = ', istep0 + write(nu_diag,1021) ' npt = ', npt + write(nu_diag,1021) ' diagfreq = ', diagfreq + write(nu_diag,1011) ' print_global = ', print_global + write(nu_diag,1011) ' print_points = ', print_points + write(nu_diag,1031) ' bfbflag = ', trim(bfbflag) + write(nu_diag,1021) ' numin = ', numin + write(nu_diag,1021) ' numax = ', numax + write(nu_diag,1033) ' histfreq = ', histfreq(:) + write(nu_diag,1023) ' histfreq_n = ', histfreq_n(:) + write(nu_diag,1011) ' hist_avg = ', hist_avg + if (.not. hist_avg) write(nu_diag,1031) ' History data will be snapshots' + write(nu_diag,1031) ' history_dir = ', trim(history_dir) + write(nu_diag,1031) ' history_file = ', trim(history_file) + write(nu_diag,1021) ' history_precision= ', history_precision + write(nu_diag,1031) ' history_format = ', trim(history_format) if (write_ic) then - write(nu_diag,*) ' Initial condition will be written in ', & + write(nu_diag,1031) ' Initial condition will be written in ', & trim(incond_dir) endif - write(nu_diag,1030) ' dumpfreq = ', & - trim(dumpfreq) - write(nu_diag,1020) ' dumpfreq_n = ', dumpfreq_n - write(nu_diag,1010) ' dump_last = ', dump_last - write(nu_diag,1010) ' restart = ', restart - write(nu_diag,*) ' restart_dir = ', & - trim(restart_dir) - write(nu_diag,*) ' restart_ext = ', restart_ext - write(nu_diag,*) ' restart_coszen = ', restart_coszen - write(nu_diag,*) ' restart_format = ', & - trim(restart_format) - write(nu_diag,*) ' lcdf64 = ', & - lcdf64 - write(nu_diag,*) ' restart_file = ', & - trim(restart_file) - write(nu_diag,*) ' pointer_file = ', & - trim(pointer_file) - write(nu_diag,*) ' use_restart_time = ', use_restart_time - write(nu_diag,*) ' ice_ic = ', & - trim(ice_ic) + write(nu_diag,1031) ' dumpfreq = ', trim(dumpfreq) + write(nu_diag,1021) ' dumpfreq_n = ', dumpfreq_n + write(nu_diag,1011) ' dump_last = ', dump_last + write(nu_diag,1011) ' restart = ', restart + write(nu_diag,1031) ' restart_dir = ', trim(restart_dir) + write(nu_diag,1011) ' restart_ext = ', restart_ext + write(nu_diag,1011) ' restart_coszen = ', restart_coszen + write(nu_diag,1031) ' restart_format = ', trim(restart_format) + write(nu_diag,1011) ' lcdf64 = ', lcdf64 + write(nu_diag,1031) ' restart_file = ', trim(restart_file) + write(nu_diag,1031) ' pointer_file = ', trim(pointer_file) + write(nu_diag,1011) ' use_restart_time = ', use_restart_time + write(nu_diag,1031) ' ice_ic = ', trim(ice_ic) if (trim(grid_type) /= 'rectangular' .or. & trim(grid_type) /= 'column') then - write(nu_diag,*) ' grid_file = ', & - trim(grid_file) - write(nu_diag,*) ' gridcpl_file = ', & - trim(gridcpl_file) - write(nu_diag,*) ' bathymetry_file = ', & - trim(bathymetry_file) - write(nu_diag,*) ' kmt_file = ', & - trim(kmt_file) - endif - write(nu_diag,1010) ' close_boundaries = ', & - close_boundaries - write(nu_diag,1010) ' orca_halogrid = ', & - orca_halogrid - - if (kdyn == 3) then - write(nu_diag,1020) ' maxits_nonlin = ', maxits_nonlin - write(nu_diag,1030) ' precond = ', precond - write(nu_diag,1020) ' dim_fgmres = ', dim_fgmres - write(nu_diag,1020) ' dim_pgmres = ', dim_pgmres - write(nu_diag,1020) ' maxits_fgmres = ', maxits_fgmres - write(nu_diag,1020) ' maxits_pgmres = ', maxits_pgmres - write(nu_diag,1010) ' monitor_nonlin = ', monitor_nonlin - write(nu_diag,1010) ' monitor_fgmres = ', monitor_fgmres - write(nu_diag,1010) ' monitor_pgmres = ', monitor_pgmres - write(nu_diag,1030) ' ortho_type = ', ortho_type - write(nu_diag,1008) ' reltol_nonlin = ', reltol_nonlin - write(nu_diag,1008) ' reltol_fgmres = ', reltol_fgmres - write(nu_diag,1008) ' reltol_pgmres = ', reltol_pgmres - write(nu_diag,1030) ' algo_nonlin = ', algo_nonlin - write(nu_diag,1010) ' use_mean_vrel = ', use_mean_vrel - if (algo_nonlin == 'anderson') then - write(nu_diag,1020) ' fpfunc_andacc = ', fpfunc_andacc - write(nu_diag,1020) ' dim_andacc = ', dim_andacc - write(nu_diag,1008) ' reltol_andacc = ', reltol_andacc - write(nu_diag,1005) ' damping_andacc = ', damping_andacc - write(nu_diag,1020) ' start_andacc = ', start_andacc - endif + write(nu_diag,1031) ' grid_file = ', trim(grid_file) + write(nu_diag,1031) ' gridcpl_file = ', trim(gridcpl_file) + write(nu_diag,1031) ' bathymetry_file = ', trim(bathymetry_file) + write(nu_diag,1031) ' kmt_file = ', trim(kmt_file) endif + write(nu_diag,1011) ' orca_halogrid = ', orca_halogrid - write(nu_diag,1010) ' conserv_check = ', conserv_check + write(nu_diag,1011) ' conserv_check = ', conserv_check - write(nu_diag,1020) ' fyear_init = ', & - fyear_init - write(nu_diag,1020) ' ycycle = ', ycycle - write(nu_diag,*) ' atm_data_type = ', & - trim(atm_data_type) + write(nu_diag,1021) ' fyear_init = ', fyear_init + write(nu_diag,1021) ' ycycle = ', ycycle + write(nu_diag,1031) ' atm_data_type = ', trim(atm_data_type) if (trim(atm_data_type) /= 'default') then - write(nu_diag,*) ' atm_data_dir = ', & - trim(atm_data_dir) - write(nu_diag,*) ' precip_units = ', & - trim(precip_units) + write(nu_diag,1031) ' atm_data_dir = ', trim(atm_data_dir) + write(nu_diag,1031) ' precip_units = ', trim(precip_units) elseif (trim(atm_data_type)=='default') then - write(nu_diag,*) ' default_season = ', trim(default_season) + write(nu_diag,1031) ' default_season = ', trim(default_season) endif if (wave_spec) then - write(nu_diag,*) ' wave_spec_file = ', trim(wave_spec_file) + write(nu_diag,1031) ' wave_spec_file = ', trim(wave_spec_file) endif if (trim(bgc_data_type) == 'ncar' .or. & trim(ocn_data_type) == 'ncar') then - write(nu_diag,*) ' oceanmixed_file = ', & - trim(oceanmixed_file) + write(nu_diag,1031) ' oceanmixed_file = ', trim(oceanmixed_file) endif if (cpl_bgc) then - write(nu_diag,1000) ' BGC coupling is switched ON' + write(nu_diag,*) 'BGC coupling is switched ON' else - write(nu_diag,1000) ' BGC coupling is switched OFF' - endif - write(nu_diag,*) ' bgc_data_type = ', & - trim(bgc_data_type) - write(nu_diag,*) ' fe_data_type = ', & - trim(fe_data_type) - write(nu_diag,*) ' ice_data_type = ', & - trim(ice_data_type) - write(nu_diag,*) ' bgc_data_dir = ', & - trim(bgc_data_dir) - write(nu_diag,*) ' ocn_data_type = ', & - trim(ocn_data_type) + write(nu_diag,*) 'BGC coupling is switched OFF' + endif + write(nu_diag,1031) ' bgc_data_type = ', trim(bgc_data_type) + write(nu_diag,1031) ' fe_data_type = ', trim(fe_data_type) + write(nu_diag,1031) ' ice_data_type = ', trim(ice_data_type) + write(nu_diag,1031) ' bgc_data_dir = ', trim(bgc_data_dir) + write(nu_diag,1031) ' ocn_data_type = ', trim(ocn_data_type) if (trim(bgc_data_type) /= 'default' .or. & trim(ocn_data_type) /= 'default') then - write(nu_diag,*) ' ocn_data_dir = ', & - trim(ocn_data_dir) - write(nu_diag,1010) ' restore_ocn = ', & - restore_ocn + write(nu_diag,1031) ' ocn_data_dir = ', trim(ocn_data_dir) + write(nu_diag,1011) ' restore_ocn = ', restore_ocn endif - write(nu_diag,1010) ' restore_ice = ', & - restore_ice + write(nu_diag,1011) ' restore_ice = ', restore_ice if (restore_ice .or. restore_ocn) & - write(nu_diag,1020) ' trestore = ', trestore + write(nu_diag,1021) ' trestore = ', trestore write(nu_diag,*) ' ' - write(nu_diag,'(a30,2f8.2)') 'Diagnostic point 1: lat, lon =', & + write(nu_diag,'(a31,2f8.2)') 'Diagnostic point 1: lat, lon =', & latpnt(1), lonpnt(1) - write(nu_diag,'(a30,2f8.2)') 'Diagnostic point 2: lat, lon =', & + write(nu_diag,'(a31,2f8.2)') 'Diagnostic point 2: lat, lon =', & latpnt(2), lonpnt(2) + write(nu_diag,*) ' ' ! tracer restarts - write(nu_diag,1010) ' restart_age = ', restart_age - write(nu_diag,1010) ' restart_FY = ', restart_FY - write(nu_diag,1010) ' restart_lvl = ', restart_lvl - write(nu_diag,1010) ' restart_pond_cesm = ', restart_pond_cesm - write(nu_diag,1010) ' restart_pond_lvl = ', restart_pond_lvl - write(nu_diag,1010) ' restart_pond_topo = ', restart_pond_topo - write(nu_diag,1010) ' restart_iso = ', restart_iso - write(nu_diag,1010) ' restart_aero = ', restart_aero - write(nu_diag,1010) ' restart_fsd = ', restart_fsd - - write(nu_diag,1020) ' n_iso = ', n_iso - write(nu_diag,1020) ' n_aero = ', n_aero - write(nu_diag,1020) ' n_zaero = ', n_zaero - write(nu_diag,1020) ' n_algae = ', n_algae - write(nu_diag,1020) ' n_doc = ', n_doc - write(nu_diag,1020) ' n_dic = ', n_dic - write(nu_diag,1020) ' n_don = ', n_don - write(nu_diag,1020) ' n_fed = ', n_fed - write(nu_diag,1020) ' n_fep = ', n_fep + write(nu_diag,1011) ' restart_age = ', restart_age + write(nu_diag,1011) ' restart_FY = ', restart_FY + write(nu_diag,1011) ' restart_lvl = ', restart_lvl + write(nu_diag,1011) ' restart_pond_cesm= ', restart_pond_cesm + write(nu_diag,1011) ' restart_pond_lvl = ', restart_pond_lvl + write(nu_diag,1011) ' restart_pond_topo= ', restart_pond_topo + write(nu_diag,1011) ' restart_iso = ', restart_iso + write(nu_diag,1011) ' restart_aero = ', restart_aero + write(nu_diag,1011) ' restart_fsd = ', restart_fsd + + write(nu_diag,1021) ' n_iso = ', n_iso + write(nu_diag,1021) ' n_aero = ', n_aero + write(nu_diag,1021) ' n_zaero = ', n_zaero + write(nu_diag,1021) ' n_algae = ', n_algae + write(nu_diag,1021) ' n_doc = ', n_doc + write(nu_diag,1021) ' n_dic = ', n_dic + write(nu_diag,1021) ' n_don = ', n_don + write(nu_diag,1021) ' n_fed = ', n_fed + write(nu_diag,1021) ' n_fep = ', n_fep + write(nu_diag,*) ' ' endif ! my_task = master_task @@ -1766,6 +1779,7 @@ subroutine input_data rfracmin_in=rfracmin, rfracmax_in=rfracmax, pndaspect_in=pndaspect, hs1_in=hs1, hp1_in=hp1, & ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, conduct_in=conduct, & a_rapid_mode_in=a_rapid_mode, Rac_rapid_mode_in=Rac_rapid_mode, & + floediam_in=floediam, hfrazilmin_in=hfrazilmin, & aspect_rapid_mode_in=aspect_rapid_mode, dSdt_slow_mode_in=dSdt_slow_mode, & phi_c_slow_mode_in=phi_c_slow_mode, phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, & wave_spec_type_in = wave_spec_type, & @@ -1785,20 +1799,17 @@ subroutine input_data if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - 1000 format (a30,2x,f9.2) ! a30 to align formatted, unformatted statements - 1002 format (a20,1x,f7.2,a) - 1005 format (a30,2x,f12.6) ! float - 1006 format (a20,2x,f10.6,a) - 1007 format (a20,2x,f6.2,a) - 1008 format (a30,2x,d13.6) ! float, exponential notation - 1009 format (a20,2x,d13.6,a) ! float, exponential notation - 1010 format (a30,2x,l6) ! logical - 1012 format (a20,2x,l3,1x,a) ! logical - 1020 format (a30,2x,i6) ! integer - 1022 format (a20,2x,i3,1x,a) ! integer - 1030 format (a30, a8) ! character - 1040 format (a30,2x,6i6) ! integer - 1050 format (a30,2x,6a6) ! character + 1000 format (a20,1x,f13.6,1x,a) ! float + 1002 format (a20,5x,f9.2,1x,a) + 1009 format (a20,1x,d13.6,1x,a) + 1010 format (a20,8x,l6,1x,a) ! logical + 1011 format (a20,1x,l6) + 1020 format (a20,8x,i6,1x,a) ! integer + 1021 format (a20,1x,i6) + 1023 format (a20,1x,6i6) + 1030 format (a20,a14,1x,a) ! character + 1031 format (a20,1x,a,a) + 1033 format (a20,1x,6a6) end subroutine input_data diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90 index 4b92c2a42..b21908e77 100644 --- a/cicecore/cicedynB/general/ice_step_mod.F90 +++ b/cicecore/cicedynB/general/ice_step_mod.F90 @@ -550,22 +550,23 @@ subroutine step_therm2 (dt, iblk) logical (kind=log_kind) :: & tr_fsd, & ! floe size distribution tracers - z_tracers + z_tracers, & ! vertical biogeochemistry + solve_zsal ! zsalinity type (block) :: & this_block ! block information for current block character(len=*), parameter :: subname = '(step_therm2)' - call icepack_query_parameters(z_tracers_out=z_tracers) + call icepack_query_parameters(z_tracers_out=z_tracers,solve_zsal_out=solve_zsal) call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - ! tcraig, nltrcr used to be the number of zbgc tracers, but it's used as a zbgc flag in icepack - if (z_tracers) then + ! nltrcr is only used as a zbgc flag in icepack (number of zbgc tracers > 0) + if (z_tracers .or. solve_zsal) then nltrcr = 1 else nltrcr = 0 @@ -851,10 +852,9 @@ subroutine step_dyn_horiz (dt) use ice_dyn_evp, only: evp use ice_dyn_eap, only: eap use ice_dyn_vp, only: implicit_solver - use ice_dyn_shared, only: kdyn, ktransport + use ice_dyn_shared, only: kdyn use ice_flux, only: init_history_dyn -!deprecate upwind use ice_transport_driver, only: advection, transport_upwind, transport_remap - use ice_transport_driver, only: advection, transport_remap + use ice_transport_driver, only: advection, transport_upwind, transport_remap real (kind=dbl_kind), intent(in) :: & dt ! dynamics time step @@ -875,13 +875,10 @@ subroutine step_dyn_horiz (dt) ! Horizontal ice transport !----------------------------------------------------------------- -!deprecate upwind if (ktransport > 0) then - if (ktransport > 0 .and. advection == 'remap') then -!deprecate upwind if (advection == 'upwind') then -!deprecate upwind call transport_upwind (dt) ! upwind -!deprecate upwind else + if (advection == 'upwind') then + call transport_upwind (dt) ! upwind + elseif (advection == 'remap') then call transport_remap (dt) ! incremental remapping -!deprecate upwind endif endif end subroutine step_dyn_horiz diff --git a/cicecore/cicedynB/infrastructure/comm/mpi/ice_communicate.F90 b/cicecore/cicedynB/infrastructure/comm/mpi/ice_communicate.F90 index a7d186083..1c369ef93 100644 --- a/cicecore/cicedynB/infrastructure/comm/mpi/ice_communicate.F90 +++ b/cicecore/cicedynB/infrastructure/comm/mpi/ice_communicate.F90 @@ -13,18 +13,6 @@ module ice_communicate use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted -#if defined key_oasis3 || key_oasis3mct - use cpl_oasis3 -#endif - -#if defined key_oasis4 - use cpl_oasis4 -#endif - -#if defined key_iomput - use lib_mpp, only: mpi_comm_opa ! MPP library -#endif - implicit none private @@ -83,13 +71,7 @@ subroutine init_communicate(mpicom) if (present(mpicom)) then ice_comm = mpicom else -#if (defined key_oasis3 || defined key_oasis3mct || defined key_oasis4) - ice_comm = localComm ! communicator from NEMO/OASISn -#elif defined key_iomput - ice_comm = mpi_comm_opa ! communicator from NEMO/XIOS -#else ice_comm = MPI_COMM_WORLD ! Global communicator -#endif endif call MPI_INITIALIZED(flag,ierr) diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 67129c911..a354efb6b 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -47,7 +47,7 @@ module ice_grid gridcpl_file , & ! input file for POP coupling grid info grid_file , & ! input file for POP grid info kmt_file , & ! input file for POP grid info - bathymetry_file, & ! input bathymetry for basalstress + bathymetry_file, & ! input bathymetry for seabed stress bathymetry_format, & ! bathymetry file format (default or pop) grid_spacing , & ! default of 30.e3m or set by user in namelist grid_type ! current options are rectangular (default), @@ -1663,7 +1663,7 @@ subroutine makemask field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) - !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block do i = 1, nx_block @@ -2343,9 +2343,9 @@ subroutine gridbox_verts(work_g,vbounds) end subroutine gridbox_verts !======================================================================= -! ocean bathymetry for grounded sea ice (basalstress) or icebergs +! ocean bathymetry for grounded sea ice (seabed stress) or icebergs ! currently hardwired for 40 levels (gx3, gx1 grids) -! should be read from a file instead (see subroutine read_basalstress_bathy) +! should be read from a file instead (see subroutine read_seabedstress_bathy) subroutine get_bathymetry @@ -2387,7 +2387,7 @@ subroutine get_bathymetry if (use_bathymetry) then - call read_basalstress_bathy + call read_seabedstress_bathy else @@ -2504,14 +2504,14 @@ end subroutine get_bathymetry_popfile !======================================================================= -! Read bathymetry data for basal stress calculation (grounding scheme for +! Read bathymetry data for seabed stress calculation (grounding scheme for ! landfast ice) in CICE stand-alone mode. When CICE is in coupled mode ! (e.g. CICE-NEMO), hwater should be uptated at each time level so that ! it varies with ocean dynamics. ! ! author: Fred Dupont, CMC - subroutine read_basalstress_bathy + subroutine read_seabedstress_bathy ! use module use ice_read_write @@ -2526,7 +2526,7 @@ subroutine read_basalstress_bathy logical (kind=log_kind) :: diag=.true. - character(len=*), parameter :: subname = '(read_basalstress_bathy)' + character(len=*), parameter :: subname = '(read_seabedstress_bathy)' if (my_task == master_task) then write (nu_diag,*) ' ' @@ -2553,7 +2553,7 @@ subroutine read_basalstress_bathy call icepack_warnings_flush(nu_diag) endif - end subroutine read_basalstress_bathy + end subroutine read_seabedstress_bathy !======================================================================= diff --git a/cicecore/drivers/direct/hadgem3/CICE.F90 b/cicecore/drivers/direct/hadgem3/CICE.F90 index 72bf1b747..e444dcd40 100644 --- a/cicecore/drivers/direct/hadgem3/CICE.F90 +++ b/cicecore/drivers/direct/hadgem3/CICE.F90 @@ -1,8 +1,8 @@ !======================================================================= -! Copyright (c) 2020, Triad National Security, LLC +! Copyright (c) 2021, Triad National Security, LLC ! All rights reserved. ! -! Copyright 2020. Triad National Security, LLC. This software was +! Copyright 2021. Triad National Security, LLC. This software was ! produced under U.S. Government contract DE-AC52-06NA25396 for Los ! Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy. The U.S. diff --git a/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 b/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 index 49cf12ce1..5f91ed584 100644 --- a/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 +++ b/cicecore/drivers/direct/hadgem3/CICE_InitMod.F90 @@ -71,7 +71,7 @@ subroutine cice_init use ice_domain, only: init_domain_blocks use ice_domain_size, only: ncat, nfsd use ice_dyn_eap, only: init_eap, alloc_dyn_eap - use ice_dyn_shared, only: kdyn, init_dyn, basalstress, alloc_dyn_shared + use ice_dyn_shared, only: kdyn, init_dyn, alloc_dyn_shared use ice_dyn_vp, only: init_vp use ice_flux, only: init_coupler_flux, init_history_therm, & init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux diff --git a/cicecore/drivers/mct/cesm1/CICE_copyright.txt b/cicecore/drivers/mct/cesm1/CICE_copyright.txt index 007cc58cf..e10da1e77 100644 --- a/cicecore/drivers/mct/cesm1/CICE_copyright.txt +++ b/cicecore/drivers/mct/cesm1/CICE_copyright.txt @@ -1,7 +1,7 @@ -! Copyright (c) 2020, Triad National Security, LLC +! Copyright (c) 2021, Triad National Security, LLC ! All rights reserved. ! -! Copyright 2020. Triad National Security, LLC. This software was +! Copyright 2021. Triad National Security, LLC. This software was ! produced under U.S. Government contract DE-AC52-06NA25396 for Los ! Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy. The U.S. diff --git a/cicecore/drivers/nuopc/cmeps/CICE_copyright.txt b/cicecore/drivers/nuopc/cmeps/CICE_copyright.txt index 007cc58cf..e10da1e77 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_copyright.txt +++ b/cicecore/drivers/nuopc/cmeps/CICE_copyright.txt @@ -1,7 +1,7 @@ -! Copyright (c) 2020, Triad National Security, LLC +! Copyright (c) 2021, Triad National Security, LLC ! All rights reserved. ! -! Copyright 2020. Triad National Security, LLC. This software was +! Copyright 2021. Triad National Security, LLC. This software was ! produced under U.S. Government contract DE-AC52-06NA25396 for Los ! Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy. The U.S. diff --git a/cicecore/drivers/nuopc/dmi/CICE.F90 b/cicecore/drivers/nuopc/dmi/CICE.F90 index ec1963d38..2fd0c9f88 100644 --- a/cicecore/drivers/nuopc/dmi/CICE.F90 +++ b/cicecore/drivers/nuopc/dmi/CICE.F90 @@ -1,8 +1,8 @@ !======================================================================= -! Copyright (c) 2020, Triad National Security, LLC +! Copyright (c) 2021, Triad National Security, LLC ! All rights reserved. ! -! Copyright 2020. Triad National Security, LLC. This software was +! Copyright 2021. Triad National Security, LLC. This software was ! produced under U.S. Government contract DE-AC52-06NA25396 for Los ! Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy. The U.S. diff --git a/cicecore/drivers/standalone/cice/CICE.F90 b/cicecore/drivers/standalone/cice/CICE.F90 index ec1963d38..2fd0c9f88 100644 --- a/cicecore/drivers/standalone/cice/CICE.F90 +++ b/cicecore/drivers/standalone/cice/CICE.F90 @@ -1,8 +1,8 @@ !======================================================================= -! Copyright (c) 2020, Triad National Security, LLC +! Copyright (c) 2021, Triad National Security, LLC ! All rights reserved. ! -! Copyright 2020. Triad National Security, LLC. This software was +! Copyright 2021. Triad National Security, LLC. This software was ! produced under U.S. Government contract DE-AC52-06NA25396 for Los ! Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy. The U.S. diff --git a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 index bd818211e..9f6f42f28 100644 --- a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 +++ b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 @@ -548,7 +548,7 @@ subroutine coupling_prep (iblk) Tref (:,:,iblk), Qref (:,:,iblk), & fresh (:,:,iblk), fsalt (:,:,iblk), & fhocn (:,:,iblk), & - fswthru (:,:,iblk), & + fswthru (:,:,iblk), & fswthru_vdr (:,:,iblk), & fswthru_vdf (:,:,iblk), & fswthru_idr (:,:,iblk), & diff --git a/cicecore/drivers/standalone/cice/CICE_RunMod.F90_debug b/cicecore/drivers/standalone/cice/CICE_RunMod.F90_debug index 8f5de17ea..5f7eebe31 100644 --- a/cicecore/drivers/standalone/cice/CICE_RunMod.F90_debug +++ b/cicecore/drivers/standalone/cice/CICE_RunMod.F90_debug @@ -247,7 +247,7 @@ plabeld = 'post step_therm2' call debug_ice (iblk, plabeld) - endif + endif ! ktherm > 0 enddo ! iblk !$OMP END PARALLEL DO @@ -394,8 +394,8 @@ albpnd, albcnt, apeff_ai, fpond, fresh, l_mpond_fresh, & alvdf_ai, alidf_ai, alvdr_ai, alidr_ai, fhocn_ai, & fresh_ai, fsalt_ai, fsalt, & - fswthru_ai, fhocn, fswthru, scale_factor, snowfrac, & - fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf, & + fswthru_ai, fhocn, scale_factor, snowfrac, & + fswthru, fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf, & swvdr, swidr, swvdf, swidf, Tf, Tair, Qa, strairxT, strairyT, & fsens, flat, fswabs, flwout, evap, Tref, Qref, & scale_fluxes, frzmlt_init, frzmlt @@ -589,11 +589,12 @@ evap (:,:,iblk), & Tref (:,:,iblk), Qref (:,:,iblk), & fresh (:,:,iblk), fsalt (:,:,iblk), & - fhocn (:,:,iblk), fswthru (:,:,iblk), & - fswthru_vdr(:,:,iblk), & - fswthru_vdf(:,:,iblk), & - fswthru_idr(:,:,iblk), & - fswthru_idf(:,:,iblk), & + fhocn (:,:,iblk), & + fswthru (:,:,iblk), & + fswthru_vdr (:,:,iblk), & + fswthru_vdf (:,:,iblk), & + fswthru_idr (:,:,iblk), & + fswthru_idf (:,:,iblk), & faero_ocn(:,:,:,iblk), & alvdr (:,:,iblk), alidr (:,:,iblk), & alvdf (:,:,iblk), alidf (:,:,iblk), & diff --git a/cicecore/version.txt b/cicecore/version.txt index 83a606cb9..e16cf8bfe 100644 --- a/cicecore/version.txt +++ b/cicecore/version.txt @@ -1 +1 @@ -CICE 6.1.3 +CICE 6.1.4 diff --git a/configuration/scripts/cice.batch.csh b/configuration/scripts/cice.batch.csh index 54a2be711..6d1f735a4 100755 --- a/configuration/scripts/cice.batch.csh +++ b/configuration/scripts/cice.batch.csh @@ -185,7 +185,15 @@ cat >> ${jobfile} << EOFB #SBATCH --qos=standby EOFB -else if (${ICE_MACHINE} =~ daley* || ${ICE_MACHINE} =~ banting*) then +else if (${ICE_MACHINE} =~ daley* || ${ICE_MACHINE} =~ banting* ) then +cat >> ${jobfile} << EOFB +#PBS -N ${ICE_CASENAME} +#PBS -j oe +#PBS -l select=${nnodes}:ncpus=${corespernode}:mpiprocs=${taskpernodelimit}:ompthreads=${nthrds} +#PBS -l walltime=${batchtime} +EOFB + +else if (${ICE_MACHINE} =~ freya* ) then cat >> ${jobfile} << EOFB #PBS -N ${ICE_CASENAME} #PBS -j oe diff --git a/configuration/scripts/cice.launch.csh b/configuration/scripts/cice.launch.csh index dfbffd6ab..a05b3a9d3 100755 --- a/configuration/scripts/cice.launch.csh +++ b/configuration/scripts/cice.launch.csh @@ -141,6 +141,18 @@ aprun -n ${ntasks} -N ${taskpernodelimit} -d ${nthrds} ./cice >&! \$ICE_RUNLOG_F EOFR endif +#======= +else if (${ICE_MACHINE} =~ freya*) then +if (${ICE_COMMDIR} =~ serial*) then +cat >> ${jobfile} << EOFR +aprun -n 1 -N 1 -d 1 ./cice >&! \$ICE_RUNLOG_FILE +EOFR +else +cat >> ${jobfile} << EOFR +aprun -n ${ntasks} -N ${taskpernodelimit} -d ${nthrds} ./cice >&! \$ICE_RUNLOG_FILE +EOFR +endif + #======= else if (${ICE_MACHINE} =~ hera*) then cat >> ${jobfile} << EOFR diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 3139726f5..f34db14f0 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -63,7 +63,7 @@ nfsd = 1 nilyr = 7 nslyr = 1 - nblyr = 7 + nblyr = 1 orca_halogrid = .false. / @@ -108,9 +108,8 @@ dSdt_slow_mode = -5.0e-8 phi_c_slow_mode = 0.05 phi_i_mushy = 0.85 - sw_redist = .false. - sw_frac = 0.9d0 - sw_dtemp = 0.02d0 + hfrazilmin = 0.05d0 + floediam = 300.0d0 / &dynamics_nml @@ -130,8 +129,9 @@ Cf = 17. Ktens = 0. e_ratio = 2. - basalstress = .false. - k1 = 8. + seabed_stress = .false. + seabed_stress_method = 'LKD' + k1 = 7.5 k2 = 15. alphab = 20. threshold_hw = 30. @@ -170,6 +170,9 @@ dT_mlt = 1.5 rsnw_mlt = 1500. kalg = 0.6 + sw_redist = .false. + sw_frac = 0.9d0 + sw_dtemp = 0.02d0 / &ponds_nml @@ -193,7 +196,7 @@ natmiter = 5 atmiter_conv = 0.0d0 ustar_min = 0.0005 - emissivity = 0.95 + emissivity = 0.985 fbot_xfer_type = 'constant' update_ocn_f = .false. l_mpond_fresh = .false. diff --git a/configuration/scripts/machines/Macros.freya_gnu b/configuration/scripts/machines/Macros.freya_gnu new file mode 100644 index 000000000..dafa4fd9b --- /dev/null +++ b/configuration/scripts/machines/Macros.freya_gnu @@ -0,0 +1,40 @@ +#============================================================================== +# Makefile macros for DMI freya +#============================================================================== +# For use with GNU compiler +#============================================================================== + +#INCLDIR := -I. -I/usr/include +#SLIBS := + +#--- Compiler/preprocessor --- +FC := ftn +CC := cc +CXX := CC +CPP := /usr/bin/cpp +CPPFLAGS := -P -traditional # ALLOW fortran double backslash "\\" +SCC := gcc +SFC := ftn + +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 +#-xHost + +FREEFLAGS := -ffree-form +FFLAGS := -fconvert=big-endian -fbacktrace -ffree-line-length-none +#-xHost + +ifeq ($(ICE_BLDDEBUG), true) + FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow +else + FFLAGS += -O2 #FROM BANTING + #FFLAGS := -O2 -ffloat-store -march=native -ffree-line-length-non # DMI BUILD +endif +LD:= $(FC) +LDFLAGS := $(FFLAGS) -v + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -fopenmp + CFLAGS += -fopenmp + FFLAGS += -fopenmp +endif diff --git a/configuration/scripts/machines/Macros.freya_intel b/configuration/scripts/machines/Macros.freya_intel new file mode 100644 index 000000000..f40ca4e23 --- /dev/null +++ b/configuration/scripts/machines/Macros.freya_intel @@ -0,0 +1,69 @@ +#============================================================================== +# Makefile macros for DMI Freya based on ECCC banting +#============================================================================== +# For use with intel compiler +#============================================================================== + +#INCLDIR := -I. -I/usr/include +#SLIBS := + +#--- Compiler/preprocessor --- +FC := ftn +CC := cc +CXX := CC +CPP := /usr/bin/cpp +CPPFLAGS := -P -traditional # ALLOW fortran double backslash "\\" +SCC := gcc +SFC := ftn + +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise +# Additional flags +FIXEDFLAGS := -132 +FREEFLAGS := -FR +FFLAGS := -fp-model source -convert big_endian -assume byterecl -ftz -traceback -no-wrap-margin +#-xHost + +ifeq ($(ICE_BLDDEBUG), true) + FFLAGS += -O0 -g -check -fpe0 -ftrapuv -fp-model except -check noarg_temp_created +# -heap-arrays 1024 +else + FFLAGS += -O2 +endif +LD := $(FC) +LDFLAGS := $(FFLAGS) -v +#ifeq ($(ICE_BLDDEBUG), true) +#FFLAGS := -O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created +#FFLAGS := -g -O0 -traceback -fp-model precise -fp-stack-check -fpe0 +#else +#FFLAGS := -r8 -i4 -O2 -align all -w -ftz -assume byterecl +# FFLAGS := -O2 -fp-model precise -assume byterecl -ftz -traceback -xHost +#endif +# Preprocessor flags +#CPPDEFS := -DLINUX $(ICE_CPPDEFS) + +# Linker flags + +# Additional flags + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp +endif + +#--- NetCDF --- +#ifeq ($(IO_TYPE), netcdf) +# +#endif +# +#ifeq ($(IO_TYPE), netcdf_bin) +# CPPDEFS := $(CPPDEFS) -Dncdf +#endif + +### if using parallel I/O, load all 3 libraries. PIO must be first! +#ifeq ($(ICE_IOTYPE), pio) +# PIO_PATH:=/usr/projects/climate/SHARED_CLIMATE/software/conejo/pio/1.7.2/intel-13.0.1/openmpi-1.6.3/netcdf-3.6.3-parallel-netcdf-1.3.1/include +# INCLDIR += -I$(PIO_PATH) +# SLIBS := $(SLIBS) -L$(PIO_PATH) -lpio +#endif diff --git a/configuration/scripts/machines/env.freya_gnu b/configuration/scripts/machines/env.freya_gnu new file mode 100755 index 000000000..b655d6dd0 --- /dev/null +++ b/configuration/scripts/machines/env.freya_gnu @@ -0,0 +1,39 @@ +#!/bin/csh -f + +set inp = "undefined" +if ($#argv == 1) then + set inp = $1 +endif + +if ("$inp" != "-nomodules") then + +source /opt/modules/default/init/csh # Initialize modules for csh + Clear environment +module rm PrgEnv-intel +module rm PrgEnv-cray +module rm PrgEnv-gnu +module add PrgEnv-gnu +#module load PrgEnv-intel # Intel compiler +#module load cray-mpich # MPI (Cray MPICH) +module add cray-netcdf # NetCDF +module add cray-hdf5 # HDF5 +#setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem +setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME freya +setenv ICE_MACHINE_MACHINFO "Cray XC50, GNU Xeon Gold 6148 (Skylake) NOT SURE-TILL" +setenv ICE_MACHINE_ENVNAME gnu +setenv ICE_MACHINE_ENVINFO "gcc/7.2.0, cray-mpich/7.7.0, cray-netcdf/4.4.1.1.6" +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR /data/${USER}/cice_original/run/ +setenv ICE_MACHINE_INPUTDATA /data/${USER}/cice_original/ +setenv ICE_MACHINE_BASELINE /data/${USER}/cice_original/dbaselines/ +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 36 # tasks per node +#setenv ICE_MACHINE_MAXRUNLENGTH 9 +setenv ICE_MACHINE_ACCT P0000000 +setenv ICE_MACHINE_QUEUE "development" +setenv ICE_MACHINE_BLDTHRDS 18 +setenv ICE_MACHINE_QSTAT "qstat " diff --git a/configuration/scripts/machines/env.freya_intel b/configuration/scripts/machines/env.freya_intel new file mode 100755 index 000000000..dcbc1f8ba --- /dev/null +++ b/configuration/scripts/machines/env.freya_intel @@ -0,0 +1,38 @@ +#!/bin/csh -f + +set inp = "undefined" +if ($#argv == 1) then + set inp = $1 +endif + +if ("$inp" != "-nomodules") then + +source /opt/modules/default/init/csh # Initialize modules for csh +# Clear environment +module rm PrgEnv-intel +module rm PrgEnv-cray +module rm PrgEnv-gnu +module add PrgEnv-intel +#module load PrgEnv-intel # Intel compiler +#module load cray-mpich # MPI (Cray MPICH) +module add cray-netcdf # NetCDF +module add cray-hdf5 # HDF5 +#setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME freya +setenv ICE_MACHINE_MACHINFO "Cray XC50, Intel Xeon Gold 6148 (Skylake) NOT SURE-TILL" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO "Intel 18.0.0.128, cray-mpich/7.7.0, cray-netcdf/4.4.1.1.6" +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR /data/${USER}/cice_original/run/ +setenv ICE_MACHINE_INPUTDATA /data/${USER}/cice_original/ +setenv ICE_MACHINE_BASELINE /data/${USER}/cice_original/dbaselines/ +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 36 # tasks per node +#setenv ICE_MACHINE_MAXRUNLENGTH 9 +setenv ICE_MACHINE_ACCT P0000000 +setenv ICE_MACHINE_QUEUE "development" +setenv ICE_MACHINE_BLDTHRDS 18 +setenv ICE_MACHINE_QSTAT "qstat " diff --git a/configuration/scripts/options/set_env.box2001 b/configuration/scripts/options/set_env.box2001 old mode 100755 new mode 100644 diff --git a/configuration/scripts/options/set_nml.alt01 b/configuration/scripts/options/set_nml.alt01 index 8b66913c6..98124b3f2 100644 --- a/configuration/scripts/options/set_nml.alt01 +++ b/configuration/scripts/options/set_nml.alt01 @@ -1,4 +1,5 @@ nilyr = 1 +use_leap_years = .false. ice_ic = 'default' restart = .false. distribution_type = 'roundrobin' @@ -15,7 +16,8 @@ kitd = 0 ktherm = 0 conduct = 'bubbly' kdyn = 0 -basalstress = .true. +seabed_stress = .true. +seabed_stress_method = 'probabilistic' use_bathymetry = .true. shortwave = 'ccsm3' albedo_type = 'constant' diff --git a/configuration/scripts/options/set_nml.alt03 b/configuration/scripts/options/set_nml.alt03 index f82491d9d..507f56a1b 100644 --- a/configuration/scripts/options/set_nml.alt03 +++ b/configuration/scripts/options/set_nml.alt03 @@ -21,5 +21,5 @@ tfrz_option = 'linear_salt' revised_evp = .false. Ktens = 0. e_ratio = 2. -basalstress = .true. +seabed_stress = .true. use_bathymetry = .true. diff --git a/configuration/scripts/options/set_nml.alt04 b/configuration/scripts/options/set_nml.alt04 index e3689fe82..937704294 100644 --- a/configuration/scripts/options/set_nml.alt04 +++ b/configuration/scripts/options/set_nml.alt04 @@ -22,7 +22,7 @@ kevp_kernel = 102 fbot_xfer_type = 'Cdn_ocn' shortwave = 'dEdd' formdrag = .true. -advection = 'remap' +advection = 'upwind' kstrength = 0 krdg_partic = 0 krdg_redist = 0 diff --git a/configuration/scripts/options/set_nml.box2001 b/configuration/scripts/options/set_nml.box2001 old mode 100755 new mode 100644 diff --git a/configuration/scripts/options/set_nml.boxrestore b/configuration/scripts/options/set_nml.boxrestore index 6789e1ff8..d00ec41c8 100644 --- a/configuration/scripts/options/set_nml.boxrestore +++ b/configuration/scripts/options/set_nml.boxrestore @@ -24,5 +24,5 @@ revised_evp = .true. kstrength = 0 krdg_partic = 0 krdg_redist = 0 -basalstress = .true. +seabed_stress = .true. restore_ice = .true. diff --git a/configuration/scripts/options/set_nml.gbox80 b/configuration/scripts/options/set_nml.gbox80 old mode 100755 new mode 100644 diff --git a/configuration/scripts/options/set_nml.gx1 b/configuration/scripts/options/set_nml.gx1 index df6e2cd2b..e1d18dc8b 100644 --- a/configuration/scripts/options/set_nml.gx1 +++ b/configuration/scripts/options/set_nml.gx1 @@ -1,5 +1,8 @@ dt = 3600.0 runtype = 'initial' +year_init = 2005 +use_leap_years = .true. +use_restart_time = .false. ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' grid_format = 'bin' grid_type = 'displaced_pole' @@ -10,11 +13,9 @@ maskhalo_dyn = .true. maskhalo_remap = .true. maskhalo_bound = .true. 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' +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx1' +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' +precip_units = 'mks' ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/COREII' bgc_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/WOA/MONTHLY' - diff --git a/configuration/scripts/options/set_nml.gx1coreii b/configuration/scripts/options/set_nml.gx1coreii new file mode 100644 index 000000000..44b334194 --- /dev/null +++ b/configuration/scripts/options/set_nml.gx1coreii @@ -0,0 +1,9 @@ +year_init = 1997 +use_leap_years = .false. +use_restart_time = .true. +fyear_init = 2005 +ycycle = 1 +atm_data_format = 'bin' +atm_data_type = 'LYq' +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/COREII' +precip_units = 'mm_per_sec' diff --git a/configuration/scripts/options/set_nml.gx3 b/configuration/scripts/options/set_nml.gx3 index 8cdb7e40a..1a2fe62a5 100644 --- a/configuration/scripts/options/set_nml.gx3 +++ b/configuration/scripts/options/set_nml.gx3 @@ -1,15 +1,18 @@ dt = 3600.0 runtype = 'initial' +year_init = 2005 +use_leap_years = .true. +use_restart_time = .false. ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc' grid_format = 'bin' grid_type = 'displaced_pole' grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin' kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin' bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc' -fyear_init = 1997 -atm_data_format = 'bin' -atm_data_type = 'ncar' -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/NCAR_bulk' +fyear_init = 2005 +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx3' +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/JRA55' +precip_units = 'mks' ocn_data_format = 'bin' ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/' - diff --git a/configuration/scripts/options/set_nml.gx3ncarbulk b/configuration/scripts/options/set_nml.gx3ncarbulk new file mode 100644 index 000000000..fbe0f7ae7 --- /dev/null +++ b/configuration/scripts/options/set_nml.gx3ncarbulk @@ -0,0 +1,9 @@ +year_init = 1997 +use_leap_years = .false. +use_restart_time = .true. +fyear_init = 1997 +atm_data_format = 'bin' +atm_data_type = 'ncar' +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/NCAR_bulk' +precip_units = 'mm_per_month' + diff --git a/configuration/scripts/options/set_nml.jra55 b/configuration/scripts/options/set_nml.jra55 deleted file mode 100755 index b35e99a6d..000000000 --- a/configuration/scripts/options/set_nml.jra55 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2005 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2005 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55_gx1' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_2008 b/configuration/scripts/options/set_nml.jra55_2008 deleted file mode 100755 index 042431fc0..000000000 --- a/configuration/scripts/options/set_nml.jra55_2008 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2008 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2008 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx1 b/configuration/scripts/options/set_nml.jra55_gx1 deleted file mode 100755 index b35e99a6d..000000000 --- a/configuration/scripts/options/set_nml.jra55_gx1 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2005 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2005 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55_gx1' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx1_2008 b/configuration/scripts/options/set_nml.jra55_gx1_2008 deleted file mode 100755 index bf92e93d1..000000000 --- a/configuration/scripts/options/set_nml.jra55_gx1_2008 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2008 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2008 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55_gx1' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx3 b/configuration/scripts/options/set_nml.jra55_gx3 deleted file mode 100755 index c2d1eefdc..000000000 --- a/configuration/scripts/options/set_nml.jra55_gx3 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2005 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2005 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55_gx3' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx3_2008 b/configuration/scripts/options/set_nml.jra55_gx3_2008 deleted file mode 100755 index 89c7dc55d..000000000 --- a/configuration/scripts/options/set_nml.jra55_gx3_2008 +++ /dev/null @@ -1,17 +0,0 @@ -year_init = 2008 -ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc' -grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin' -kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin' -bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc' -use_leap_years = .true. -use_restart_time = .false. -maskhalo_dyn = .true. -maskhalo_remap = .true. -maskhalo_bound = .true. -fyear_init = 2008 -atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/JRA55' -atm_data_format = 'nc' -atm_data_type = 'JRA55_gx3' -precip_units = 'mks' -ocn_data_dir = 'default' -bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.kevp102 b/configuration/scripts/options/set_nml.kevp102 new file mode 100644 index 000000000..3a5dc3dbd --- /dev/null +++ b/configuration/scripts/options/set_nml.kevp102 @@ -0,0 +1 @@ +kevp_kernel = 102 diff --git a/configuration/scripts/options/set_nml.yi2008 b/configuration/scripts/options/set_nml.yi2008 new file mode 100644 index 000000000..ab7a0ecee --- /dev/null +++ b/configuration/scripts/options/set_nml.yi2008 @@ -0,0 +1,3 @@ +year_init = 2008 +fyear_init = 2008 + diff --git a/configuration/scripts/options/set_nml.zsal b/configuration/scripts/options/set_nml.zsal new file mode 100644 index 000000000..c099fed8c --- /dev/null +++ b/configuration/scripts/options/set_nml.zsal @@ -0,0 +1,9 @@ +nblyr = 7 +ktherm = 1 +sw_redist = .true. +tfrz_option = 'linear_salt' +tr_brine = .true. +solve_zsal = .true. +ice_ic = 'default' +restart = .false. + diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 386c29e41..1ed489730 100755 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -5,9 +5,9 @@ smoke gx3 1x4 debug,diag1,run2day smoke gx3 4x1 debug,diag1,run5day restart gx3 8x2 debug smoke gx3 8x2 diag24,run1year,medium -decomp gx3 4x2x25x29x5 -smoke gx3 4x2 diag1,run5day smoke_gx3_8x2_diag1_run5day -smoke gx3 4x1 diag1,run5day,thread smoke_gx3_8x2_diag1_run5day +decomp gx3 4x2x25x29x5 none +smoke gx3 4x2 diag1,run5day smoke_gx3_8x2_diag1_run5day +smoke gx3 4x1 diag1,run5day,thread smoke_gx3_8x2_diag1_run5day restart gx1 40x4 droundrobin,medium restart tx1 40x4 dsectrobin,medium restart gx3 4x4 none @@ -33,21 +33,26 @@ smoke gbox80 1x1 boxslotcyl smoke gx3 8x2 bgcz smoke gx3 8x2 bgcz,debug smoke gx3 8x1 bgcskl,debug -#smoke gx3 4x1 bgcz,thread smoke_gx3_8x2_bgcz +#smoke gx3 4x1 bgcz,thread smoke_gx3_8x2_bgcz restart gx1 4x2 bgcsklclim,medium restart gx1 8x1 bgczclim,medium -smoke gx1 24x1 jra55_gx1_2008,medium,run90day -smoke gx3 8x1 jra55_gx3_2008,medium,run90day -restart gx1 24x1 jra55_gx1,short -restart gx3 8x1 jra55_gx3,short +smoke gx1 24x1 medium,run90day,yi2008 +smoke gx3 8x1 medium,run90day,yi2008 +restart gx1 24x1 short +restart gx3 8x1 short smoke gx3 4x2 fsd1,diag24,run5day,debug smoke gx3 8x2 fsd12,diag24,run5day,short restart gx3 4x2 fsd12,debug,short smoke gx3 8x2 fsd12ww3,diag24,run1day,medium smoke gx3 4x1 isotope,debug restart gx3 8x2 isotope -restart gx3 4x4 iobinary +restart gx3 4x4 gx3ncarbulk,iobinary restart gx3 4x4 histall,precision8,cdf64 smoke gx3 30x1 bgcz,histall smoke gx3 14x2 fsd12,histall smoke gx3 4x1 dynpicard,medium +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 diff --git a/configuration/scripts/tests/decomp_suite.ts b/configuration/scripts/tests/decomp_suite.ts index 1e1ae3112..4eb5394d9 100644 --- a/configuration/scripts/tests/decomp_suite.ts +++ b/configuration/scripts/tests/decomp_suite.ts @@ -1,6 +1,6 @@ # Test Grid PEs Sets BFB-compare restart gx3 4x2x25x29x4 dslenderX2 -decomp gx3 4x2x25x29x5 +decomp gx3 4x2x25x29x5 none sleep 30 restart gx3 1x1x50x58x4 droundrobin,thread restart_gx3_4x2x25x29x4_dslenderX2 restart gx3 4x1x25x116x1 dslenderX1,thread restart_gx3_4x2x25x29x4_dslenderX2 diff --git a/configuration/scripts/tests/first_suite.ts b/configuration/scripts/tests/first_suite.ts index b06ee7e0b..bf6c813f6 100644 --- a/configuration/scripts/tests/first_suite.ts +++ b/configuration/scripts/tests/first_suite.ts @@ -1,5 +1,5 @@ -# Test Grid PEs Sets BFB-compare -smoke gx3 8x2 diag1,run5day +# Test Grid PEs Sets BFB-compare +smoke gx3 8x2 diag1,run5day restart gx3 4x2x25x29x4 dslenderX2 logbfb gx3 4x2x25x29x4 dslenderX2,diag1,reprosum -smoke gx3 1x2 run2day +smoke gx3 1x2 run2day diff --git a/configuration/scripts/tests/io_suite.ts b/configuration/scripts/tests/io_suite.ts index 3e98642e9..a17e3f625 100755 --- a/configuration/scripts/tests/io_suite.ts +++ b/configuration/scripts/tests/io_suite.ts @@ -1,15 +1,16 @@ # Test Grid PEs Sets BFB-compare # some iobinary configurations fail due to bathymetry netcdf file requirement, remove them -restart gx3 8x4 debug,histall,iobinary,precision8 -#restart gx3 12x2 alt01,histall,iobinary -restart gx3 16x2 alt02,histall,iobinary,precision8 -#restart gx3 4x2 alt03,histall,iobinary -restart gx3 8x4 alt04,histall,iobinary,precision8 -restart gx3 4x4 alt05,histall,iobinary -restart gx3 32x1 bgcz,histall,iobinary,precision8 -restart gx3 16x2 bgcskl,histall,iobinary -restart gx3 14x2 isotope,histall,iobinary,precision8 -restart gx3 16x2 fsd12,histall,iobinary +# iobinary cannot work with JRA55 because netcdf is turned off +restart gx3 8x4 gx3ncarbulk,debug,histall,iobinary,precision8 +#restart gx3 12x2 gx3ncarbulk,alt01,histall,iobinary +restart gx3 16x2 gx3ncarbulk,alt02,histall,iobinary,precision8 +#restart gx3 4x2 gx3ncarbulk,alt03,histall,iobinary +restart gx3 8x4 gx3ncarbulk,alt04,histall,iobinary,precision8 +restart gx3 4x4 gx3ncarbulk,alt05,histall,iobinary +restart gx3 32x1 gx3ncarbulk,bgcz,histall,iobinary,precision8 +restart gx3 16x2 gx3ncarbulk,bgcskl,histall,iobinary +restart gx3 14x2 gx3ncarbulk,isotope,histall,iobinary,precision8 +restart gx3 16x2 gx3ncarbulk,fsd12,histall,iobinary restart gx3 32x1 debug,histall,ionetcdf restart gx3 15x2 alt01,histall,ionetcdf,precision8,cdf64 diff --git a/configuration/scripts/tests/nothread_quicksuite.ts b/configuration/scripts/tests/nothread_quicksuite.ts index 535347438..997b59b28 100644 --- a/configuration/scripts/tests/nothread_quicksuite.ts +++ b/configuration/scripts/tests/nothread_quicksuite.ts @@ -1,5 +1,4 @@ # Test Grid PEs Sets BFB-compare - restart gx3 16x1 diag1 smoke gx3 1x1 debug,diag1,run2day smoke gx3 4x1 debug,diag1,run2day,thread diff --git a/configuration/scripts/tests/nothread_suite.ts b/configuration/scripts/tests/nothread_suite.ts index 49f834a98..afe1963b3 100644 --- a/configuration/scripts/tests/nothread_suite.ts +++ b/configuration/scripts/tests/nothread_suite.ts @@ -14,7 +14,7 @@ smoke gx3 16x1 diag24,run1year,medium #restart tx1 160x1 dsectrobin,medium restart gx3 16x1 none -restart gx3 16x1 iobinary +restart gx3 16x1 gx3ncarbulk,iobinary restart gx3 12x1 alt01 restart gx3 16x1 alt02 @@ -43,11 +43,11 @@ smoke gbox128 24x1 boxrestore,short,debug restart gbox80 1x1 box2001 smoke gbox80 1x1 boxslotcyl -smoke gx3 16x1 jra55_gx3_2008,medium,run90day -restart gx3 12x1 jra55_gx3,short +smoke gx3 16x1 medium,run90day,yi2008 +restart gx3 12x1 short #tcraig, hangs nodes intermittently on izumi -#smoke gx1 24x1 jra55_gx1_2008,medium,run90day -#restart gx1 24x1 jra55_gx1,short +#smoke gx1 24x1 medium,run90day,yi2008 +#restart gx1 24x1 short smoke gx3 16x1 bgcz smoke gx3 16x1 bgcz,debug @@ -56,7 +56,7 @@ smoke gx3 24x1 bgcskl,debug #restart gx1 128x1 bgcsklclim,medium #restart gx1 256x1 bgczclim,medium -decomp gx3 8x1x5x29x20 +decomp gx3 8x1x5x29x20 none restart gx3 1x1x50x58x4 droundrobin restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 4x1x25x116x1 dslenderX1 restart_gx3_8x1x25x29x2_dslenderX2 restart gx3 12x1x4x29x9 dspacecurve restart_gx3_8x1x25x29x2_dslenderX2 diff --git a/configuration/scripts/tests/quick_suite.ts b/configuration/scripts/tests/quick_suite.ts index 10595e9ad..9384f0333 100644 --- a/configuration/scripts/tests/quick_suite.ts +++ b/configuration/scripts/tests/quick_suite.ts @@ -1,6 +1,6 @@ # Test Grid PEs Sets BFB-compare -smoke gx3 8x2 diag1,run5day -smoke gx3 1x1 diag1,run1day +smoke gx3 8x2 diag1,run5day +smoke gx3 1x1 diag1,run1day restart gbox128 8x1 diag1 restart gx3 4x2 debug,diag1,run5day smoke gx3 4x1 diag1,run5day,thread smoke_gx3_8x2_diag1_run5day diff --git a/configuration/scripts/tests/reprosum_suite.ts b/configuration/scripts/tests/reprosum_suite.ts index 34cf51a80..d65370e0a 100644 --- a/configuration/scripts/tests/reprosum_suite.ts +++ b/configuration/scripts/tests/reprosum_suite.ts @@ -4,8 +4,8 @@ logbfb gx3 4x2x25x29x4 dslenderX2,diag1,reprosum 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 -logbfb gx3 8x2x8x10x20 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum -logbfb gx3 6x2x50x58x1 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum +logbfb gx3 1x20x5x29x80 dsectrobin,diag1,short,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum +logbfb gx3 8x2x8x10x20 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum +logbfb gx3 6x2x50x58x1 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum logbfb gx3 6x2x4x29x18 dspacecurve,diag1,maskhalo,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum -#logbfb gx3 8x2x8x10x20 droundrobin,diag1 logbfb_gx3_4x2x25x29x4_diag1_dslenderX2 +#logbfb gx3 8x2x8x10x20 droundrobin,diag1 logbfb_gx3_4x2x25x29x4_diag1_dslenderX2 diff --git a/configuration/scripts/tests/test_decomp.script b/configuration/scripts/tests/test_decomp.script index 70d7cd24a..40b84d08e 100644 --- a/configuration/scripts/tests/test_decomp.script +++ b/configuration/scripts/tests/test_decomp.script @@ -75,10 +75,16 @@ foreach decomp (${decomps}) set bfbstatus = $status if ($bfbstatus == 0) then set grade = PASS - echo "bfb baseline and test dataset are identical" + echo "bfbcomp baseline and test dataset are identical" + else if ( ${bfbstatus} == 1 ) then + set grade = FAIL + echo "bfbcomp baseline and test dataset are different" + else if ( ${bfbstatus} == 2 ) then + set grade = FAIL + echo "bfbcomp baseline missing data" else set grade = FAIL - echo "bfbcomp and test dataset are different" + echo "bfbcomp script failed" endif echo "$grade ${ICE_TESTNAME}_${decomp} bfbcomp ${base_case}" >> ${ICE_CASEDIR}/test_output endif diff --git a/configuration/scripts/tests/travis_suite.ts b/configuration/scripts/tests/travis_suite.ts index 6f475f4c5..b452bae0f 100644 --- a/configuration/scripts/tests/travis_suite.ts +++ b/configuration/scripts/tests/travis_suite.ts @@ -3,5 +3,5 @@ smoke gx3 1x2 run2day smoke gx3 1x1 debug,run1day smoke gx3 2x2 debug,run1day smoke gx3 2x1 run2day,thread smoke_gx3_1x2_run2day -restart gx3 2x1 -restart gx3 1x2 +restart gx3 2x1 +restart gx3 1x2 diff --git a/configuration/scripts/timeseries.csh b/configuration/scripts/timeseries.csh index cdd025efc..b6b3fcf2e 100755 --- a/configuration/scripts/timeseries.csh +++ b/configuration/scripts/timeseries.csh @@ -13,9 +13,9 @@ endif set basename = `echo $1 | sed -e 's#/$##' | sed -e 's/^\.\///'` # Set x-axis limits - # Manuallyl set x-axis limits + # Manually set x-axis limits #set xrange = 'set xrange ["19980101":"19981231"]' - # Let gnuplot determine x-alis limits + # Let gnuplot determine x-axis limits set xrange = '' # Determine if BASELINE dataset exists @@ -59,8 +59,10 @@ endif # Loop through each field and create the plot foreach field ($fieldlist:q) + # Add backslashes before (, ), and ^ for grep searches + set search_name = "`echo '$field' | sed 's/(/\\(/' | sed 's/)/\\)/' | sed 's/\^/\\^/'`" set fieldname = `echo "$field" | sed -e 's/([^()]*)//g'` - set search = "'$fieldname'\|istep1" + set search = "'$search_name'\|istep1" rm -f data.txt foreach line ("`egrep $search $logfile`") if ("$line" =~ *"istep1"*) then diff --git a/configuration/scripts/timeseries.py b/configuration/scripts/timeseries.py index 2b50c373a..2c36cea73 100755 --- a/configuration/scripts/timeseries.py +++ b/configuration/scripts/timeseries.py @@ -53,7 +53,7 @@ def get_data(logfile,field): # Build the regular expression to extract the data field_regex = field.replace('(','\(').replace('^','\^').replace(')','\)') number_regex = '[-+]?\d+\.?\d+([eE][-+]?\d+)?' - my_regex = '{}\s+=\s+({})\s+({})'.format(field_regex,number_regex,number_regex) + my_regex = '^{}\s+=\s+({})\s+({})'.format(field_regex,number_regex,number_regex) dtg = [] arctic = [] @@ -95,9 +95,10 @@ def plot_timeseries(log, field, dtg, arctic, antarctic, expon, dtg_base=None, ar Plot the timeseries data from the CICE log file ''' - casename = os.path.abspath(log).rstrip('/').rstrip('/logs').split('/')[-1] + import re + casename = re.sub(r"/logs", "", os.path.abspath(log).rstrip('/')).split('/')[-1] if base_dir: - base_casename = os.path.abspath(base_dir).rstrip('/').rstrip('/logs').split('/')[-1] + base_casename = re.sub(r"/logs", "", os.path.abspath(base_dir).rstrip('/')).split('/')[-1] # Load the plotting libraries, but set the logging level for matplotlib # to WARNING so that matplotlib debugging info is not printed when running diff --git a/doc/requirements.txt b/doc/requirements.txt index 8788d6ac3..218e06d1f 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -1,5 +1,5 @@ # # -sphinxcontrib-bibtex +sphinxcontrib-bibtex<2.0.0 # # diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 8ea16261d..59ddc4122 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -188,7 +188,7 @@ either Celsius or Kelvin units). "e11, e12, e22", "strain rate tensor components", "" "ecci", "yield curve minor/major axis ratio, squared", "1/4" "eice(n)", "energy of melting of ice per unit area (in category n)", "J/m\ :math:`^2`" - "emissivity", "emissivity of snow and ice", "0.95" + "emissivity", "emissivity of snow and ice", "0.985" "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`" @@ -553,6 +553,8 @@ either Celsius or Kelvin units). "s11, s12, s22", "stress tensor components", "" "saltmax", "max salinity, at ice base (:cite:`Bitz99`)", "3.2 ppt" "scale_factor", "scaling factor for shortwave radiation components", "" + "seabed_stress", "if true, calculate seabed stress", "F" + "seabed_stress_method", "method for calculating seabed stress (‘LKD’ or ‘probabilistic’)", "LKD" "sec", "seconds elasped into idate", "" "secday", "number of seconds in a day", "86400." "shcoef", "transfer coefficient for sensible heat", "" diff --git a/doc/source/conf.py b/doc/source/conf.py index 8d0df9777..e876980ab 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -54,7 +54,7 @@ # General information about the project. project = u'CICE' -copyright = u'2020, Triad National Security, LLC (code) and National Center for Atmospheric Research (documentation)' +copyright = u'2021, Triad National Security, LLC (code) and National Center for Atmospheric Research (documentation)' author = u'CICE-Consortium' # The version info for the project you're documenting, acts as replacement for @@ -62,9 +62,9 @@ # built documents. # # The short X.Y version. -version = u'6.1.3' +version = u'6.1.4' # The full version, including alpha/beta/rc tags. -version = u'6.1.3' +version = u'6.1.4' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/developer_guide/dg_dynamics.rst b/doc/source/developer_guide/dg_dynamics.rst index eac19b1f6..c94d47b35 100644 --- a/doc/source/developer_guide/dg_dynamics.rst +++ b/doc/source/developer_guide/dg_dynamics.rst @@ -50,9 +50,9 @@ abort if set. To override the abort, use value 102 for testing. Transport ----------------- -The transport (advection) methods are found in **cicecore/cicedynB/dynamics/**. Only the incremental -remapping method is supported at this time, and is set in namelist via the ``advection`` variable. -Transport can be turned off by setting ``advection = none`` or ``ktransport = -1``. +The transport (advection) methods are found in **cicecore/cicedynB/dynamics/**. Two methods are supported, +upwind and remap. These are set in namelist via the ``advection`` variable. +Transport can be disabled with the ``ktransport`` namelist variable. Infrastructure diff --git a/doc/source/developer_guide/dg_forcing.rst b/doc/source/developer_guide/dg_forcing.rst new file mode 100644 index 000000000..90ef843b0 --- /dev/null +++ b/doc/source/developer_guide/dg_forcing.rst @@ -0,0 +1,217 @@ +:tocdepth: 3 + +.. _forcing: + +Standalone Forcing +====================== + +Users are strongly encouraged to run CICE in a coupled system (see :ref:`coupl`) to improve +quality of science. The standalone mode is best used for technical testing +and only preliminary science testing. Several different input forcing datasets have +been implemented over the history of CICE. Some have become obsolete, others +have been supplanted by newer forcing data options, and others have been implemented +by outside users and are not testable by the Consortium. The forcing code has +generally not been maintained by the Consortium and only a subset of the code is +tested by the Consortium. + +The forcing implementation can be found in the file +**cicecore/cicedynB/general/ice_forcing.F90**. As noted above, only a subset of the +forcing modes are tested and supported. In many ways, the implemetation is fairly +primitive, in part due to historical reasons and in part because standalone runs +are discouraged for evaluating complex science. In general, most implementations +use aspects of the following approach, + +- Input files are organized by year. +- Namelist inputs ``fyear`` and ``ycycle`` specify the forcing year dataset. +- The forcing year is computed on the fly and is assumed to be cyclical over the forcing dataset length defined by ``ycycle``. +- The namelist ``atm_dat_dir`` specifies the directory of the atmosphere input data files and the namelist ``atm_data_type`` defines the atmospheric forcing mode. +- The namelist ``ocn_dat_dir`` specifies the directory of the ocean input data files and the namelist ``ocn_data_type`` defines the ocean forcing mode. +- The filenames follow a particular naming convention that is defined in the source code (ie. subroutine **JRA55_gx1_files**). The forcing year is typically found just before the **.nc** part of the filename and there are tools (subroutine **file_year**) to update the filename based on the model year and appropriate forcing year. +- The input data time axis is generally NOT read by the forcing subroutine. The forcing frequency is hardwired into the model and the file record number is computed based on the forcing frequency and model time. Mixing leap year input data and noleap model calendars (and vice versa) is not handled particularly gracefully. The CICE model does not read or check against the input data time axis. +- Data is read on the model grid, no spatial interpolation exists. +- Data is often time interpolated linearly between two input timestamps to the model time each model timestep. + +In general, the following variables need to be defined by the forcing module, + +From Atmosphere: + +- zlvl = atmosphere level height (m) +- uatm = model grid i-direction wind velocity component (m/s) +- vatm = model grid j-direction wind velocity component (m/s) +- strax = model grid i-direction wind stress (N/m^2) +- stray = model grid j-direction wind stress (N/m^2) +- potT = air potential temperature (K) +- Tair = air temperature (K) +- Qa = specific humidity (kg/kg) +- rhoa = air density (kg/m^3) +- flw = incoming longwave radiation (W/m^2) +- fsw = incoming shortwave radiation (W/m^2) +- swvdr = sw down, visible, direct (W/m^2) +- swvdf = sw down, visible, diffuse (W/m^2) +- swidr = sw down, near IR, direct (W/m^2) +- swidf = sw down, near IR, diffuse (W/m^2) +- frain = rainfall rate (kg/m^2 s) +- fsnow = snowfall rate (kg/m^2 s) + +From Ocean: + +- uocn = ocean current, x-direction (m/s) +- vocn = ocean current, y-direction (m/s) +- ss_tltx = sea surface slope, x-direction (m/m) +- ss_tlty = sea surface slope, y-direction (m/m) +- hwater = water depth for basal stress calc (landfast ice) +- sss = sea surface salinity (ppt) +- sst = sea surface temperature (C) +- frzmlt = freezing/melting potential (W/m^2) +- frzmlt_init= frzmlt used in current time step (W/m^2) +- Tf = freezing temperature (C) +- qdp = deep ocean heat flux (W/m^2), negative upward +- hmix = mixed layer depth (m) +- daice_da= data assimilation concentration increment rate (concentration s-1)(only used in hadgem drivers) + +All variables have reasonable but static defaults and these will be used in ``default`` mode. + +To advance the forcing, the subroutines **get_forcing_atmo** and +**get_forcing_ocn** are called each timestep from the step +loop. That subroutine computes the forcing year (``fyear``), calls the appropriate +forcing data method, and then calls **prepare_forcing** which converts the +input data fields to model forcing fields. + +.. _JRA55forcing: + +JRA55 Atmosphere Forcing +------------------------- + +The current default atmosphere forcing for gx3, gx1, and tx1 standalone grids for +Consortium testing is the JRA55 forcing +dataset :cite:`Tsujino18`. The Consortium has released 5 years of forcing data, +2005-2009 for gx3, gx1, and tx1 grids. Each year is a separate file and +the dataset is on a gregorian time axis which includes leap days. + +.. _fig-jra55data: + +.. figure:: ./figures/jra55data.png + :align: center + :scale: 100% + + Schematic of JRA55 CICE forcing file generation. + +The forcing is read and interpolated in subroutine **JRA55_data**. In particular, +air temperature (``airtmp``), east and north wind speed (``wndewd`` and ``wndnwd``), +specific humidity (``spchmd``), incoming short and longwave radiation (``glbrad`` and ``dswsfc``), +and precipitation (``ttlpcp``) are read from the input files. The JRA55 reanalysis is +run with updated initial conditions every 6 hours and output is written every 3 hours. +The four state fields (air temperature, winds, specific humidity) +are instantaneous data, while the three flux fields (radition, precipitation) are 3 +hour averages. In the JRA55 forcing files provided by the Consortium, the time +defined for 3 hour average fields is shifted 3 hours to the start time of the 3 +hour interval. **NOTE that this is different +from the implementation on the original JRA55 files and also different from how models +normally define time on an accumulated/averaged field**. This is all shown +schematically in Figure :ref:`fig-jra55data`. + +The state fields are linearly +time interpolated between input timestamps +while the flux fields are read and held constant during each 3 hour model period. +The forcing frequency is hardwired to 3 hours in the implementation, +and the record number is computed based on the time of the current model year. +Time interpolation coefficients are computed in the **JRA55_data** subroutine. + +The forcing data is converted to model inputs in the subroutine **prepare_forcing** +called in **get_forcing_atmo**. To clarify, the JRA55 input data includes + +- uatm = model grid i-direction wind velocity component (m/s) +- vatm = model grid j-direction wind velocity component (m/s) +- Tair = air temperature (K) +- Qa = specific humidity (kg/kg) +- flw = incoming longwave radiation (W/m^2) +- fsw = incoming shortwave radiation (W/m^2) +- fsnow = snowfall rate (kg/m^2 s) + +and model forcing inputs are derived from those fields and the defaults. + +Because the input files are on the gregorian time axis, the model can run with the regular +365 day (noleap) calendar, but in that case, the Feb 29 input data will be used on +March 1, and all data +after March 1 will be shifted one day. December 31 in leap years will be skipped when +running with a CICE calendar with no leap days. + + +.. _NCARforcing: + +NCAR Atmosphere Forcing +------------------------- + +The NCAR atmospheric forcing was used in earlier standalone runs on the gx3 grid, and the +Consortium continues to do some limited testing with this forcing dataset. +Monthly average data for fsw, cldf, fsnow are read. 6-hourly data for +Tair, uatm, vatm, rhoa, and Qa are also read. +Users are encouraged to switch to the JRA55 (see :ref:`JRA55forcing`) dataset. This +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 +---------------------------- + +The default atmosphere forcing option sets the atmosphere forcing +internally. No files are read. Values for forcing fields are defined +at initialization in subroutine **init_coupler_flux** and held +constant thereafter. Different conditions can be specified thru the +``default_season`` namelist variable. + + +.. _box2001forcing: + +Box2001 Atmosphere Forcing +------------------------- + +The box2001 forcing dataset in generated internally. No files are read. The +dataset is used to test an idealized box case as defined in :cite:`Hunke01`. + + +.. _otheratmforcing: + +Other Atmosphere Forcing +------------------------- + +There are a few other atmospheric forcing modes, as defined by ``atm_data_type``, but +they are not tested by the Consortium on a regular basis. + + +.. _defaultocnforcing: + +Default Ocean Forcing +------------------------- + +The ``default`` ocean setting is the standard setting used in standalone CICE runs. +In this mode, the sea surface salinity is set to 34 ppt and the sea surface +temperature is set to the freezing temperature at all grid points and +held constant unless the mixed layer parameterization is turned on, in which +case the SST evolves. Other ocean coupling fields are set to zero. No files are read. + + +.. _otherocnforcing: + +Other Ocean Forcing +------------------------- + +There are a few other ocean forcing modes, as defined by ``ocn_data_type``, but +they are not tested by the Consortium on a regular basis. + diff --git a/doc/source/developer_guide/figures/jra55data.png b/doc/source/developer_guide/figures/jra55data.png new file mode 100644 index 000000000..5fee7b062 Binary files /dev/null and b/doc/source/developer_guide/figures/jra55data.png differ diff --git a/doc/source/developer_guide/index.rst b/doc/source/developer_guide/index.rst index 1b016f6d9..ab5b2d1e6 100644 --- a/doc/source/developer_guide/index.rst +++ b/doc/source/developer_guide/index.rst @@ -14,6 +14,7 @@ Developer Guide dg_about.rst dg_dynamics.rst dg_driver.rst + dg_forcing.rst dg_icepack.rst dg_scripts.rst dg_other.rst diff --git a/doc/source/intro/copyright.rst b/doc/source/intro/copyright.rst index 671140b3f..f09f6c58d 100644 --- a/doc/source/intro/copyright.rst +++ b/doc/source/intro/copyright.rst @@ -5,7 +5,7 @@ Copyright ============================= -© Copyright 2020, Triad National Security LLC. All rights reserved. +© Copyright 2021, Triad National Security LLC. All rights reserved. This software was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department diff --git a/doc/source/intro/quickstart.rst b/doc/source/intro/quickstart.rst index 547d6ef20..56d19ee70 100644 --- a/doc/source/intro/quickstart.rst +++ b/doc/source/intro/quickstart.rst @@ -16,6 +16,9 @@ You will probably have to download some inputdata, see the `CICE wiki h_{cu}`. The maximum seabed stress depends on the weight of the ridge above hydrostatic balance and the value of :math:`k_2`. It is, however, the parameter :math:`k_1` that has the most notable impact on the simulated extent of landfast ice. -The value of :math:`k_1` can be changed at runtime using the namelist variable ``k1``. The grounding scheme can be turned on or off using the namelist logical basalstress. +The value of :math:`k_1` can be changed at runtime using the namelist variable ``k1``. -Note that the user must provide a bathymetry field for using this grounding -scheme. It is suggested to have a bathymetry field with water depths larger than -5 m that represents well shallow water regions such as the Laptev Sea and the -East Siberian Sea. To prevent unrealistic grounding, :math:`T_b` is set to zero when :math:`h_{wu}` +To prevent unrealistic grounding, :math:`T_b` is set to zero when :math:`h_{wu}` is larger than 30 m. This maximum value is chosen based on observations of large keels in the Arctic Ocean :cite:`Amundrud04`. - +Seabed stress based on probabilistic approach +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This more sophisticated grounding parameterization computes the seabed stress based +on the probability of contact between the ice thickness distribution +(ITD) and the seabed. Multi-thickness category models such as CICE typically use a +few thickness categories (5-10). This crude representation of the ITD +does not resolve the tail of the ITD, which is crucial for grounding +events. + +To represent the tail of the distribution, the simulated ITD is +converted to a positively skewed probability function :math:`f(x)` +with :math:`x` the sea ice thickness. The mean and variance are set +equal to the ones of the original ITD. A +log-normal distribution is used for :math:`f(x)`. + +It is assumed that the bathymetry :math:`y` (at the 't' point) follows a normal +distribution :math:`b(y)`. The mean of :math:`b(y)` comes from the user's bathymetry field and the +standard deviation :math:`\sigma_b` is currently fixed to 2.5 m. Two +possible improvements would be to specify a distribution based on high +resolution bathymetry data and to take into account variations of the +water depth due to changes in the sea surface height. + +Assuming hydrostatic balance and neglecting the impact of snow, the draft of floating ice of thickness +:math:`x` is :math:`D(x)=\rho_i x / \rho_w` where :math:`\rho_i` is the sea ice density. Hence, the probability of contact (:math:`P_c`) between the +ITD and the seabed is given by + +.. math:: + P_c=\int_{0}^{\inf} \int_{0}^{D(x)} g(x)b(y) dy dx \label{prob_contact}. + +:math:`T_b` is first calculated at the 't' point (referred to as :math:`T_{bt}`). :math:`T_{bt}` depends on the weight of the ridge in excess of hydrostatic balance. The parameterization first calculates + +.. math:: + T_{bt}^*=\mu_s g \int_{0}^{\inf} \int_{0}^{D(x)} (\rho_i x - \rho_w + y)g(x)b(y) dy dx, \\ + :label: Tbt + +and then obtains :math:`T_{bt}` by multiplying :math:`T_{bt}^*` by :math:`e^{-\alpha_b * (1 - a_i)}` (similar to what is done for ``seabed_stress_method`` = ``LKD``). + +To calculate :math:`T_{bt}^*` in equation :eq:`Tbt`, :math:`f(x)` and :math:`b(y)` are discretized using many small categories (100). :math:`f(x)` is discretized between 0 and 50 m while :math:`b(y)` is truncated at plus and minus three :math:`\sigma_b`. :math:`f(x)` is also modified by setting it to zero after a certain percentile of the log-normal distribution. This percentile, which is currently set to 99.7%, notably affects the simulation of landfast ice and is used as a tuning parameter. Its impact is similar to the one of the parameter :math:`k_1` for the LKD method. + +:math:`T_b` at the 'u' point is calculated from the 't' point values around it according to + +.. math:: + T_b=\max[T_{bt}(i,j),T_{bt}(i+1,j),T_{bt}(i,j+1),T_{bt}(i+1,j+1)]. \\ + :label: Tb + +Following again the LKD method, the seabed stress coefficients are finally expressed as + +.. math:: + C_b= T_b (\sqrt{u^2+v^2}+u_0)^{-1}, \\ + :label: Cb2 + .. _internal-stress: *************** @@ -311,7 +382,11 @@ 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``). +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: @@ -342,9 +417,7 @@ 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. -The ice strength :math:`P` -is a function of the ice thickness and concentration -as described in the `Icepack Documentation `_. 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``. +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``. .. _stress-evp: diff --git a/doc/source/science_guide/sg_horiztrans.rst b/doc/source/science_guide/sg_horiztrans.rst index 33b37564e..bafb4c72f 100644 --- a/doc/source/science_guide/sg_horiztrans.rst +++ b/doc/source/science_guide/sg_horiztrans.rst @@ -33,7 +33,7 @@ introductory comments in **ice\_transport\_remap.F90**. Prognostic equations for ice and/or snow density may be included in future model versions but have not yet been implemented. -One transport scheme is available, the incremental +Two transport schemes are available: upwind and the incremental remapping scheme of :cite:`Dukowicz00` as modified for sea ice by :cite:`Lipscomb04`. The remapping scheme has several desirable features: diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 227a63663..ccf7f0356 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -33,10 +33,6 @@ can be found in :ref:`cicecpps`. The following CPPs are available. "ESMF_INTERFACE", "Turns on ESMF support in a subset of driver code. Also USE_ESMF_LIB and USE_ESMF_METADATA" "FORTRANUNDERSCORE", "Used in ice_shr_reprosum86.c to support Fortran-C interfaces. This should generally be turned on at all times. There are other CPPs (FORTRANDOUBULEUNDERSCORE, FORTRANCAPS, etc) in ice_shr_reprosum.c that are generally not used in CICE but could be useful if problems arise in the Fortran-C interfaces" "GPTL", "Turns on GPTL initialization if needed for PIO" - "key_oasis3", "Leverages Oasis CPPs to define the local MPI communicator" - "key_oasis3mct", "Leverages Oasis CPPs to define the local MPI communicator" - "key_oasis4", "Leverages Oasis CPPs to define the local MPI communicator" - "key_iomput", "Leverages Oasis CPPs to define the local MPI communicator" "NO_F2003", "Turns off some Fortran 2003 features" "NO_I8", "Converts integer*8 to integer*4. This could have adverse affects for certain algorithms including the ddpdd implementation associated with the ``bfbflag``" "NO_R16", "Converts real*16 to real*8. This could have adverse affects for certain algorithms including the lsum16 implementation associated with the ``bfbflag``" @@ -216,7 +212,6 @@ grid_nml "``bathymetry_file``", "string", "name of bathymetry file to be read", "‘unknown_bathymetry_file’" "``bathymetry_format``", "``default``", "NetCDF depth field", "‘default’" "", "``pop``", "pop thickness file in cm in ascii format", "" - "``close_boundaries``", "logical", "set land on edges of grid", "``.false.``" "``dxrect``", "real", "x-direction grid spacing for rectangular grid in cm", "0.0" "``dyrect``", "real", "y-direction grid spacing for rectangular grid in cm", "0.0" "``gridcpl_file``", "string", "input file for coupling grid info", "'unknown_gridcpl_file'" @@ -239,7 +234,7 @@ grid_nml "``nilyr``", "integer", "number of vertical layers in ice", "0" "``nslyr``", "integer", "number of vertical layers in snow", "0" "``orca_halogrid``", "logical", "use orca haloed grid for data/grid read", "``.false.``" - "``use_bathymetry``", "logical", "use read in bathymetry file for basalstress option", "``.false.``" + "``use_bathymetry``", "logical", "use read in bathymetry file for seabedstress option", "``.false.``" "", "", "", "" domain_nml @@ -334,6 +329,8 @@ thermo_nml "``conduct``", "``bubbly``", "conductivity scheme :cite:`Pringle07`", "``bubbly``" "", "``MU71``", "conductivity :cite:`Maykut71`", "" "``dSdt_slow_mode``", "real", "slow drainage strength parameter m/s/K", "-1.5e-7" + "``floediam``", "real", "effective floe diameter for lateral melt in m", "300.0" + "``hfrazilmin``", "real", "min thickness of new frazil ice in m", "0.05" "``kitd``", "``0``", "delta function ITD approximation", "1" "", "``1``", "linear remapping ITD approximation", "" "``ksno``", "real", "snow thermal conductivity", "0.3" @@ -344,9 +341,6 @@ thermo_nml "``phi_c_slow_mode``", ":math:`0<\phi_c < 1`", "critical liquid fraction", "0.05" "``phi_i_mushy``", ":math:`0<\phi_i < 1`", "solid fraction at lower boundary", "0.85" "``Rac_rapid_mode``", "real", "critical Rayleigh number", "10.0" - "``sw_redist``", "logical", "redistribute internal shortwave to surface", "``.false.``" - "``sw_frac``", "real", "fraction redistributed", "0.9" - "``sw_dtemp``", "real", "temperature difference from melt to start redistributing", "0.02" "", "", "", "" .. _dynamics_nml: @@ -360,11 +354,10 @@ dynamics_nml "", "", "", "" "``advection``", "``remap``", "linear remapping advection scheme", "``remap``" - "", "``none``", "advection off", "" + "", "``upwind``", "donor cell advection", "" "``alphab``", "real", ":math:`\alpha_{b}` factor in :cite:`Lemieux16`", "20.0" "``arlx``", "real", "revised_evp value", "300.0" "``brlx``", "real", "revised_evp value", "300.0" - "``basalstress``", "logical", "use basal stress parameterization for landfast ice", "``.false.``" "``Cf``", "real", "ratio of ridging work to PE change in ridging", "17.0" "``coriolis``", "``constant``", "constant coriolis value = 1.46e-4 s\ :math:`^{-1}`", "``latitude``" "", "``latitude``", "coriolis variable by latitude", "" @@ -411,6 +404,9 @@ dynamics_nml "``reltol_fgmres``", "real", "relative tolerance for FGMRES solver", "1e-2" "``reltol_pgmres``", "real", "relative tolerance for PGMRES preconditioner", "1e-6" "``revised_evp``", "logical", "use revised EVP formulation", "``.false.``" + "``seabed_stress``", "logical", "use seabed stress parameterization for landfast ice", "``.false.``" + "``seabed_stress_method``", "``LKD``", "linear keel draft method :cite:`Lemieux16`", "``LKD``" + "", "``probabilistic``", "probability of contact method (Dupont et al., in prep)", "" "``ssh_stress``", "``coupled``", "computed from coupled sea surface height gradient", "``geostrophic``" "", "``geostropic``", "computed from ocean velocity", "" "``threshold_hw``", "real", "Max water depth for grounding (see :cite:`Amundrud04`)", "30." @@ -440,6 +436,9 @@ shortwave_nml "``R_snw``", "real", "tuning parameter for snow (broadband albedo) from Delta-Eddington shortwave", "1.5" "``shortwave``", "``ccsm3``", "NCAR CCSM3 shortwave distribution method", "``ccsm3``" "", "``dEdd``", "Delta-Eddington method", "" + "``sw_dtemp``", "real", "temperature difference from melt to start redistributing", "0.02" + "``sw_frac``", "real", "fraction redistributed", "0.9" + "``sw_redist``", "logical", "redistribute internal shortwave to surface", "``.false.``" "", "", "", "" ponds_nml @@ -495,7 +494,7 @@ forcing_nml "``calc_Tsfc``", "logical", "calculate surface temperature", "``.true.``" "``default_season``", "``summer``", "forcing initial summer values", "``winter``" "", "``winter``", "forcing initial winter values", "" - "``emissivity``", "real", "emissivity of snow and ice", "0.95" + "``emissivity``", "real", "emissivity of snow and ice", "0.985" "``fbot_xfer_type``", "``Cdn_ocn``", "variabler ocean heat transfer coefficient scheme", "``constant``" "", "``constant``", "constant ocean heat transfer coefficient", "" "``fe_data_type``", "``clim``", "ocean climatology forcing value for iron", "``default``" diff --git a/doc/source/user_guide/ug_running.rst b/doc/source/user_guide/ug_running.rst index 957cfc4fc..541fa81a4 100644 --- a/doc/source/user_guide/ug_running.rst +++ b/doc/source/user_guide/ug_running.rst @@ -143,11 +143,8 @@ Some hints: - To change namelist, manually edit the **ice_in** file - To change batch settings, manually edit the top of the **cice.run** or **cice.test** (if running a test) file -- When the run scripts are submitted, the current **ice_in**, **cice.settings**, and **env.[machine]** files are copied from the case directory into the run directory. Users should generally not edit files in the run directory as these are overwritten when following the standard workflow. **cice.settings** can be sourced to establish the case values in the login shell. An alias like the following can be established to quickly switch between case and run directories:: - - alias cdrun 'cd `\grep "setenv ICE_RUNDIR" cice.settings | awk "{print "\$"NF}"`' - alias cdcase 'cd `\grep "setenv ICE_CASEDIR" cice.settings | awk "{print "\$"NF}"`' - +- When the run scripts are submitted, the current **ice_in**, **cice.settings**, and **env.[machine]** files are copied from the case directory into the run directory. Users should generally not edit files in the run directory as these are overwritten when following the standard workflow. **cice.settings** can be sourced to establish the case values in the login shell. +- Some useful aliases can be found in the :ref:`aliases` section - To turn on the debug compiler flags, set ``ICE_BLDDEBUG`` in **cice.setttings** to true. It is also possible to use the ``debug`` option (``-s debug``) when creating the case with **cice.setup** to set this option automatically. - To change compiler options, manually edit the Macros file. To add user defined preprocessor macros, modify ``ICE_CPPDEFS`` in **cice.settings** using the syntax ``-DCICE_MACRO``. - To clean the build before each compile, set ``ICE_CLEANBUILD`` in **cice.settings** to true (this is the default value), or use the ``buildclean`` option (``-s buildclean``) when creating the case with **cice.setup**. To not clean before the build, set ``ICE_CLEANBUILD`` in **cice.settings** to false, or use the ``buildincremental`` option (``-s buildincremental``) when creating the case with **cice.setup**. It is recommended that the ``ICE_CLEANBUILD`` be set to true if there are any questions about whether the build is proceeding properly. @@ -496,8 +493,9 @@ in the **env.[machine]** file. This can also be manually changed in the **cice. .. _laptops: -Porting to Laptop or Personal Computers +Porting to Laptops or Personal Computers ----------------------------------------- + To get the required software necessary to build and run CICE, and use the plotting and quality control scripts included in the repository, a `conda `_ environment file is available at : ``configuration/scripts/machines/environment.yml``. @@ -837,6 +835,50 @@ modify the scripts and input settings in the case directory, NOT the run directo In general, files in the run directory are overwritten by versions in the case directory when the model is built, submitted, and run. +.. _aliases: + +Use of Shell Aliases +------------------------- + +This section provides a list of some potentially useful shell aliases that leverage the CICE +scripts. These are not defined by CICE and are not required for using CICE. They +are provided as an example of what can be done by users. +The current **ice_in**, **cice.settings**, and **env.[machine]** files are copied from +the case directory into the run directory when the model is run. Users can create aliases +leveraging the variables in these files. Aliases like the following can be established +in shell startup files or otherwise at users discretion: + +.. code-block:: bash + + #!/bin/tcsh + # From a case or run directory, source the necessary environment files to run CICE + alias cice_env 'source env.*; source cice.settings' + # Go from case directory to run directory and back (see https://stackoverflow.com/a/34874698/) + alias cdrun 'set rundir=`\grep "setenv ICE_RUNDIR" cice.settings | awk "{print "\$"NF}"` && cd $rundir' + alias cdcase 'set casedir=`\grep "setenv ICE_CASEDIR" cice.settings | awk "{print "\$"NF}"` && cd $casedir' + + #!/bin/bash + # From case/test directory, go to run directory + alias cdrun='cd $(cice_var ICE_RUNDIR)' + # From run directory, go to case/test directory + alias cdcase='cd $(cice_var ICE_CASEDIR)' + # monitor current cice run (from ICE_RUNDIR directory) + alias cice_tail='tail -f $(ls -1t cice.runlog.* |head -1)' + # open log from last CICE run (from ICE_CASEDIR directory) + alias cice_lastrun='$EDITOR $(ls -1t logs/cice.runlog.* |head -1)' + # open log from last CICE build (from ICE_CASEDIR directory) + alias cice_lastbuild='$EDITOR $(ls -1t logs/cice.bldlog.* |head -1)' + # show CICE run directory when run in the case directory + alias cice_rundir='cice_var ICE_RUNDIR' + # open a tcsh shell and source env.* and cice.settings (useful for launching CICE in a debugger) + alias cice_shell='tcsh -c "cice_env; tcsh"' + + ## Functions + # Print the value of a CICE variable ($1) from cice.settings + cice_var() { + \grep "setenv $1" cice.settings | awk "{print "\$"3}" + } + .. _timeseries: Timeseries Plotting diff --git a/doc/source/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index d7e4a9fa4..61aa1c05f 100644 --- a/doc/source/user_guide/ug_testing.rst +++ b/doc/source/user_guide/ug_testing.rst @@ -1013,6 +1013,7 @@ Below is an example of a step-by-step procedure for testing a code change that m # Run a full regression test to verify bit-for-bit # Create a baseline dataset (only necessary if no baseline exists on the system) + # if you want to replace an existing baseline, you should first delete the directory cice.my.baseline in ${ICE_BASELINE}. # git clone the baseline code ./cice.setup -m onyx -e intel --suite base_suite --testid base0 --bgen cice.my.baseline @@ -1072,4 +1073,3 @@ If the regression comparisons fail, then you may want to run the QC test, INFO:__main__: INFO:__main__:Quality Control Test PASSED - diff --git a/icepack b/icepack index db2a47789..1fbea24e1 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit db2a4778970ae340b6bdd62eb03f60cd37a13f75 +Subproject commit 1fbea24e1c364a02dea7068e977b4e1355aef917