diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 000000000..069709970 --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.F90 diff=fortran diff --git a/.github/workflows/test-icepack.yml b/.github/workflows/test-icepack.yml index 841c372df..883031a00 100644 --- a/.github/workflows/test-icepack.yml +++ b/.github/workflows/test-icepack.yml @@ -17,7 +17,7 @@ on: defaults: run: - shell: /bin/csh {0} + shell: /bin/csh -e {0} jobs: build: @@ -83,9 +83,9 @@ jobs: conda env create -f configuration/scripts/machines/environment.yml - name: check conda env run: | - conda activate icepack && which mpicc && which mpifort && which make - mpifort --version - mpicc --version + conda activate icepack && which clang && which gfortran && which make + gfortran --version + clang --version make --version - name: check setup case run: | @@ -104,7 +104,7 @@ jobs: - name: download input data run: | cd $HOME/icepack-dirs/input - wget https://zenodo.org/record/3728287/files/Icepack_data-20200326.tar.gz && tar xvfz ICEPACK_data-20200326.tar.gz + wget --progress=dot:giga https://zenodo.org/record/3728287/files/Icepack_data-20200326.tar.gz && tar xvfz ICEPACK_data-20200326.tar.gz pwd ls -alR # - name: run case diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 1cc99db12..94d30c518 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -62,7 +62,7 @@ subroutine atmo_boundary_layer (sfctype, & Qa_iso, Qref_iso, & iso_flag, & uvel, vvel, & - Uref ) + Uref, zlvs ) use icepack_parameters, only: highfreq, natmiter, atmiter_conv @@ -118,6 +118,9 @@ subroutine atmo_boundary_layer (sfctype, & real (kind=dbl_kind), intent(out) :: & Uref ! reference height wind speed (m/s) + real (kind=dbl_kind), intent(in), optional :: & + zlvs ! atm level height (scalar quantities) (m) + ! local variables integer (kind=int_kind) :: & @@ -125,7 +128,6 @@ subroutine atmo_boundary_layer (sfctype, & real (kind=dbl_kind) :: & TsfK , & ! surface temperature in Kelvin (K) - xqq , & ! temporary variable psimh , & ! stability function at zlvl (momentum) tau , & ! stress at zlvl fac , & ! interpolation factor @@ -153,13 +155,15 @@ subroutine atmo_boundary_layer (sfctype, & re , & ! sqrt of exchange coefficient (water) rh , & ! sqrt of exchange coefficient (heat) vmag , & ! surface wind magnitude (m/s) - alz , & ! ln(zlvl /z10) + alzm , & ! ln(zlvl /z10) + alzs , & ! ln(zlvs /z10) (if zlvs present) thva , & ! virtual temperature (K) cp , & ! specific heat of moist air - hol , & ! H (at zlvl ) over L + holm , & ! H (at zlvl ) over L + hols , & ! H (at zlvs ) over L (if zlvs present) stable, & ! stability factor cpvir , & ! defined as cp_wv/cp_air - 1. - psixh ! stability function at zlvl (heat and water) + psixh ! stability function at zlvl (at zlvs if present) (heat and water) real (kind=dbl_kind), parameter :: & zTrf = c2 ! reference height for air temp (m) @@ -247,7 +251,12 @@ subroutine atmo_boundary_layer (sfctype, & thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K) delq = Qa - ssq ! spec hum dif (kg/kg) - alz = log(zlvl/zref) + alzm = log(zlvl/zref) + if (present(zlvs)) then + alzs = log(zlvs/zref) + else + alzs = alzm + endif cp = cp_air*(c1 + cpvir*ssq) !------------------------------------------------------------ @@ -275,28 +284,24 @@ subroutine atmo_boundary_layer (sfctype, & ustar_prev = ustar ! compute stability & evaluate all stability functions - hol = vonkar * gravit * zlvl & - * (tstar/thva & - + qstar/(c1/zvir+Qa)) & - / ustar**2 - hol = sign( min(abs(hol),c10), hol ) - stable = p5 + sign(p5 , hol) - xqq = max(sqrt(abs(c1 - c16*hol)) , c1) - xqq = sqrt(xqq) - - ! Jordan et al 1999 - psimhs = -(0.7_dbl_kind*hol & - + 0.75_dbl_kind*(hol-14.3_dbl_kind) & - * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind) - psimh = psimhs*stable & - + (c1 - stable)*psimhu(xqq) - psixh = psimhs*stable & - + (c1 - stable)*psixhu(xqq) + holm = compute_stability_parameter(zlvl , thva , & + ustar, tstar, & + qstar, Qa) + if (present(zlvs)) then + hols = compute_stability_parameter(zlvs , thva , & + ustar, tstar, & + qstar, Qa) + else + hols = holm + endif + + call compute_stability_function('momentum', holm, stable, psimh) + call compute_stability_function('scalar' , hols, stable, psixh) ! shift all coeffs to measurement height and stability - rd = rdn / (c1+rdn/vonkar*(alz-psimh)) - rh = rhn / (c1+rhn/vonkar*(alz-psixh)) - re = ren / (c1+ren/vonkar*(alz-psixh)) + rd = rdn / (c1+rdn/vonkar*(alzm-psimh)) + rh = rhn / (c1+rhn/vonkar*(alzs-psixh)) + re = ren / (c1+ren/vonkar*(alzs-psixh)) ! update ustar, tstar, qstar using updated, shifted coeffs ustar = rd * vmag @@ -361,16 +366,14 @@ subroutine atmo_boundary_layer (sfctype, & ! Compute diagnostics: 2m ref T, Q, U !------------------------------------------------------------ - hol = hol*zTrf/zlvl - xqq = max( c1, sqrt(abs(c1-c16*hol)) ) - xqq = sqrt(xqq) - psix2 = -c5*hol*stable + (c1-stable)*psixhu(xqq) + hols = hols*zTrf/zlvl + psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols) fac = (rh/vonkar) & - * (alz + al2 - psixh + psix2) + * (alzs + al2 - psixh + psix2) Tref = potT - delt*fac Tref = Tref - p01*zTrf ! pot temp to temp correction fac = (re/vonkar) & - * (alz + al2 - psixh + psix2) + * (alzs + al2 - psixh + psix2) Qref = Qa - delq*fac if (highfreq .and. sfctype(1:3)=='ice') then @@ -853,7 +856,7 @@ subroutine icepack_atm_boundary(sfctype, & Cdn_atm_ratio_n, & Qa_iso, Qref_iso, & uvel, vvel, & - Uref) + Uref, zlvs) character (len=3), intent(in) :: & sfctype ! ice or ocean @@ -864,7 +867,7 @@ subroutine icepack_atm_boundary(sfctype, & uatm , & ! x-direction wind speed (m/s) vatm , & ! y-direction wind speed (m/s) wind , & ! wind speed (m/s) - zlvl , & ! atm level height (m) + zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) @@ -893,7 +896,8 @@ subroutine icepack_atm_boundary(sfctype, & real (kind=dbl_kind), optional, intent(in) :: & uvel , & ! x-direction ice speed (m/s) - vvel ! y-direction ice speed (m/s) + vvel , & ! y-direction ice speed (m/s) + zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), optional, intent(out) :: & Uref ! reference height wind speed (m/s) @@ -965,7 +969,7 @@ subroutine icepack_atm_boundary(sfctype, & Qa_iso=l_Qa_iso, & Qref_iso=l_Qref_iso, & uvel=l_uvel, vvel=l_vvel, & - Uref=l_Uref) + Uref=l_Uref, zlvs=zlvs) if (icepack_warnings_aborted(subname)) return endif ! atmbndy @@ -981,31 +985,117 @@ subroutine icepack_atm_boundary(sfctype, & end subroutine icepack_atm_boundary +!======================================================================= + + function compute_stability_parameter(zlvl , thva , & + ustar, tstar, & + qstar, Qa) & + result(hol) + + real (kind=dbl_kind), intent(in) :: & + zlvl , & ! atm level height (m) + thva , & ! virtual temperature (K) + ustar , & ! turbulent scale for momentum + tstar , & ! turbulent scale for temperature + qstar , & ! turbulent scale for humidity + Qa ! specific humidity (kg/kg) + + real (kind=dbl_kind) :: & + hol ! H (at zlvl) over L + + character(len=*),parameter :: subname='(compute_stability_parameter)' + + hol = vonkar * gravit * zlvl & + * (tstar/thva & + + qstar/(c1/zvir+Qa)) & + / ustar**2 + hol = sign( min(abs(hol),c10), hol) + + end function compute_stability_parameter + +!======================================================================= + + subroutine compute_stability_function(qty, hol, stable, psi) + + character (len=*), intent(in) :: & + qty ! 'momentum' or 'scalar' + + real (kind=dbl_kind), intent(in) :: & + hol ! H over L + + real (kind=dbl_kind), intent(out) :: & + psi , & ! stability function at hol + stable ! unit step function at hol + + ! local variables + + real (kind=dbl_kind) :: & + psi_stable , & ! stable stability funcion at hol + psi_unstable ! unstable stability funcion at hol + + character(len=*),parameter :: subname='(compute_stability_function)' + + stable = p5 + sign(p5 , hol) + + psi_stable = -(0.7_dbl_kind*hol & + + 0.75_dbl_kind*(hol-14.3_dbl_kind) & + * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind) + + if(trim(qty) == 'momentum') then + psi_unstable = psi_momentum_unstable(hol) + elseif(trim(qty) == 'scalar') then + psi_unstable = psi_scalar_unstable(hol) + else + call icepack_warnings_add(subname//' incorrect qty: ' // qty) + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + endif + + psi = psi_stable*stable + (c1 - stable)*psi_unstable + + end subroutine compute_stability_function + !------------------------------------------------------------ ! Define functions !------------------------------------------------------------ !======================================================================= - real(kind=dbl_kind) function psimhu(xd) + real(kind=dbl_kind) function psi_momentum_unstable(hol) + + real(kind=dbl_kind), intent(in) :: hol + + real(kind=dbl_kind) :: xd + + xd = capital_X(hol) + + psi_momentum_unstable = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) & + - c2*atan(xd) + pih + + end function psi_momentum_unstable + +!======================================================================= + + real(kind=dbl_kind) function psi_scalar_unstable(hol) + + real(kind=dbl_kind), intent(in) :: hol + + real(kind=dbl_kind) :: xd - real(kind=dbl_kind), intent(in) :: xd + xd = capital_X(hol) - psimhu = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) & - - c2*atan(xd) + pih -!ech - c2*atan(xd) + 1.571_dbl_kind + psi_scalar_unstable = c2 * log((c1 + xd*xd)/c2) - end function psimhu + end function psi_scalar_unstable !======================================================================= - real(kind=dbl_kind) function psixhu(xd) + real(kind=dbl_kind) function capital_X(hol) - real(kind=dbl_kind), intent(in) :: xd + real(kind=dbl_kind), intent(in) :: hol - psixhu = c2 * log((c1 + xd*xd)/c2) + capital_X = sqrt(max(sqrt(abs(c1 - c16*hol)) , c1)) - end function psixhu + end function capital_X !======================================================================= diff --git a/columnphysics/icepack_orbital.F90 b/columnphysics/icepack_orbital.F90 index 2b7e9b6d6..165f3c5a1 100644 --- a/columnphysics/icepack_orbital.F90 +++ b/columnphysics/icepack_orbital.F90 @@ -212,7 +212,7 @@ end subroutine compute_coszen #ifndef CESMCOUPLED SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & - & obliqr , lambm0, mvelpp, log_print) + obliqr , lambm0, mvelpp, log_print) !------------------------------------------------------------------------------- ! @@ -261,176 +261,176 @@ SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & ! rate (arc seconds/year), phase (degrees). real (dbl_kind), parameter :: obamp(poblen) = & ! amplitudes for obliquity cos series - & (/ -2462.2214466_dbl_kind, -857.3232075_dbl_kind, -629.3231835_dbl_kind, & - & -414.2804924_dbl_kind, -311.7632587_dbl_kind, 308.9408604_dbl_kind, & - & -162.5533601_dbl_kind, -116.1077911_dbl_kind, 101.1189923_dbl_kind, & - & -67.6856209_dbl_kind, 24.9079067_dbl_kind, 22.5811241_dbl_kind, & - & -21.1648355_dbl_kind, -15.6549876_dbl_kind, 15.3936813_dbl_kind, & - & 14.6660938_dbl_kind, -11.7273029_dbl_kind, 10.2742696_dbl_kind, & - & 6.4914588_dbl_kind, 5.8539148_dbl_kind, -5.4872205_dbl_kind, & - & -5.4290191_dbl_kind, 5.1609570_dbl_kind, 5.0786314_dbl_kind, & - & -4.0735782_dbl_kind, 3.7227167_dbl_kind, 3.3971932_dbl_kind, & - & -2.8347004_dbl_kind, -2.6550721_dbl_kind, -2.5717867_dbl_kind, & - & -2.4712188_dbl_kind, 2.4625410_dbl_kind, 2.2464112_dbl_kind, & - & -2.0755511_dbl_kind, -1.9713669_dbl_kind, -1.8813061_dbl_kind, & - & -1.8468785_dbl_kind, 1.8186742_dbl_kind, 1.7601888_dbl_kind, & - & -1.5428851_dbl_kind, 1.4738838_dbl_kind, -1.4593669_dbl_kind, & - & 1.4192259_dbl_kind, -1.1818980_dbl_kind, 1.1756474_dbl_kind, & - & -1.1316126_dbl_kind, 1.0896928_dbl_kind/) + (/ -2462.2214466_dbl_kind, -857.3232075_dbl_kind, -629.3231835_dbl_kind, & + -414.2804924_dbl_kind, -311.7632587_dbl_kind, 308.9408604_dbl_kind, & + -162.5533601_dbl_kind, -116.1077911_dbl_kind, 101.1189923_dbl_kind, & + -67.6856209_dbl_kind, 24.9079067_dbl_kind, 22.5811241_dbl_kind, & + -21.1648355_dbl_kind, -15.6549876_dbl_kind, 15.3936813_dbl_kind, & + 14.6660938_dbl_kind, -11.7273029_dbl_kind, 10.2742696_dbl_kind, & + 6.4914588_dbl_kind, 5.8539148_dbl_kind, -5.4872205_dbl_kind, & + -5.4290191_dbl_kind, 5.1609570_dbl_kind, 5.0786314_dbl_kind, & + -4.0735782_dbl_kind, 3.7227167_dbl_kind, 3.3971932_dbl_kind, & + -2.8347004_dbl_kind, -2.6550721_dbl_kind, -2.5717867_dbl_kind, & + -2.4712188_dbl_kind, 2.4625410_dbl_kind, 2.2464112_dbl_kind, & + -2.0755511_dbl_kind, -1.9713669_dbl_kind, -1.8813061_dbl_kind, & + -1.8468785_dbl_kind, 1.8186742_dbl_kind, 1.7601888_dbl_kind, & + -1.5428851_dbl_kind, 1.4738838_dbl_kind, -1.4593669_dbl_kind, & + 1.4192259_dbl_kind, -1.1818980_dbl_kind, 1.1756474_dbl_kind, & + -1.1316126_dbl_kind, 1.0896928_dbl_kind/) real (dbl_kind), parameter :: obrate(poblen) = & ! rates for obliquity cosine series - & (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, & - & 31.983787_dbl_kind, 44.828336_dbl_kind, 30.973257_dbl_kind, & - & 43.668246_dbl_kind, 32.246691_dbl_kind, 30.599444_dbl_kind, & - & 42.681324_dbl_kind, 43.836462_dbl_kind, 47.439436_dbl_kind, & - & 63.219948_dbl_kind, 64.230478_dbl_kind, 1.010530_dbl_kind, & - & 7.437771_dbl_kind, 55.782177_dbl_kind, 0.373813_dbl_kind, & - & 13.218362_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, & - & 76.438310_dbl_kind, 45.815258_dbl_kind, 8.448301_dbl_kind, & - & 56.792707_dbl_kind, 49.747842_dbl_kind, 12.058272_dbl_kind, & - & 75.278220_dbl_kind, 65.241008_dbl_kind, 64.604291_dbl_kind, & - & 1.647247_dbl_kind, 7.811584_dbl_kind, 12.207832_dbl_kind, & - & 63.856665_dbl_kind, 56.155990_dbl_kind, 77.448840_dbl_kind, & - & 6.801054_dbl_kind, 62.209418_dbl_kind, 20.656133_dbl_kind, & - & 48.344406_dbl_kind, 55.145460_dbl_kind, 69.000539_dbl_kind, & - & 11.071350_dbl_kind, 74.291298_dbl_kind, 11.047742_dbl_kind, & - & 0.636717_dbl_kind, 12.844549_dbl_kind/) + (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, & + 31.983787_dbl_kind, 44.828336_dbl_kind, 30.973257_dbl_kind, & + 43.668246_dbl_kind, 32.246691_dbl_kind, 30.599444_dbl_kind, & + 42.681324_dbl_kind, 43.836462_dbl_kind, 47.439436_dbl_kind, & + 63.219948_dbl_kind, 64.230478_dbl_kind, 1.010530_dbl_kind, & + 7.437771_dbl_kind, 55.782177_dbl_kind, 0.373813_dbl_kind, & + 13.218362_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, & + 76.438310_dbl_kind, 45.815258_dbl_kind, 8.448301_dbl_kind, & + 56.792707_dbl_kind, 49.747842_dbl_kind, 12.058272_dbl_kind, & + 75.278220_dbl_kind, 65.241008_dbl_kind, 64.604291_dbl_kind, & + 1.647247_dbl_kind, 7.811584_dbl_kind, 12.207832_dbl_kind, & + 63.856665_dbl_kind, 56.155990_dbl_kind, 77.448840_dbl_kind, & + 6.801054_dbl_kind, 62.209418_dbl_kind, 20.656133_dbl_kind, & + 48.344406_dbl_kind, 55.145460_dbl_kind, 69.000539_dbl_kind, & + 11.071350_dbl_kind, 74.291298_dbl_kind, 11.047742_dbl_kind, & + 0.636717_dbl_kind, 12.844549_dbl_kind/) real (dbl_kind), parameter :: obphas(poblen) = & ! phases for obliquity cosine series - & (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, & - & 292.7252_dbl_kind, 15.3747_dbl_kind, 263.7951_dbl_kind, & - & 308.4258_dbl_kind, 240.0099_dbl_kind, 222.9725_dbl_kind, & - & 268.7809_dbl_kind, 316.7998_dbl_kind, 319.6024_dbl_kind, & - & 143.8050_dbl_kind, 172.7351_dbl_kind, 28.9300_dbl_kind, & - & 123.5968_dbl_kind, 20.2082_dbl_kind, 40.8226_dbl_kind, & - & 123.4722_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, & - & 267.2772_dbl_kind, 55.0196_dbl_kind, 152.5268_dbl_kind, & - & 49.1382_dbl_kind, 204.6609_dbl_kind, 56.5233_dbl_kind, & - & 200.3284_dbl_kind, 201.6651_dbl_kind, 213.5577_dbl_kind, & - & 17.0374_dbl_kind, 164.4194_dbl_kind, 94.5422_dbl_kind, & - & 131.9124_dbl_kind, 61.0309_dbl_kind, 296.2073_dbl_kind, & - & 135.4894_dbl_kind, 114.8750_dbl_kind, 247.0691_dbl_kind, & - & 256.6114_dbl_kind, 32.1008_dbl_kind, 143.6804_dbl_kind, & - & 16.8784_dbl_kind, 160.6835_dbl_kind, 27.5932_dbl_kind, & - & 348.1074_dbl_kind, 82.6496_dbl_kind/) + (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, & + 292.7252_dbl_kind, 15.3747_dbl_kind, 263.7951_dbl_kind, & + 308.4258_dbl_kind, 240.0099_dbl_kind, 222.9725_dbl_kind, & + 268.7809_dbl_kind, 316.7998_dbl_kind, 319.6024_dbl_kind, & + 143.8050_dbl_kind, 172.7351_dbl_kind, 28.9300_dbl_kind, & + 123.5968_dbl_kind, 20.2082_dbl_kind, 40.8226_dbl_kind, & + 123.4722_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, & + 267.2772_dbl_kind, 55.0196_dbl_kind, 152.5268_dbl_kind, & + 49.1382_dbl_kind, 204.6609_dbl_kind, 56.5233_dbl_kind, & + 200.3284_dbl_kind, 201.6651_dbl_kind, 213.5577_dbl_kind, & + 17.0374_dbl_kind, 164.4194_dbl_kind, 94.5422_dbl_kind, & + 131.9124_dbl_kind, 61.0309_dbl_kind, 296.2073_dbl_kind, & + 135.4894_dbl_kind, 114.8750_dbl_kind, 247.0691_dbl_kind, & + 256.6114_dbl_kind, 32.1008_dbl_kind, 143.6804_dbl_kind, & + 16.8784_dbl_kind, 160.6835_dbl_kind, 27.5932_dbl_kind, & + 348.1074_dbl_kind, 82.6496_dbl_kind/) ! Cosine/sine series data for computation of eccentricity and fixed vernal ! equinox longitude of perihelion (fvelp): amplitude, ! rate (arc seconds/year), phase (degrees). real (dbl_kind), parameter :: ecamp (pecclen) = & ! ampl for eccen/fvelp cos/sin series - & (/ 0.01860798_dbl_kind, 0.01627522_dbl_kind, -0.01300660_dbl_kind, & - & 0.00988829_dbl_kind, -0.00336700_dbl_kind, 0.00333077_dbl_kind, & - & -0.00235400_dbl_kind, 0.00140015_dbl_kind, 0.00100700_dbl_kind, & - & 0.00085700_dbl_kind, 0.00064990_dbl_kind, 0.00059900_dbl_kind, & - & 0.00037800_dbl_kind, -0.00033700_dbl_kind, 0.00027600_dbl_kind, & - & 0.00018200_dbl_kind, -0.00017400_dbl_kind, -0.00012400_dbl_kind, & - & 0.00001250_dbl_kind/) + (/ 0.01860798_dbl_kind, 0.01627522_dbl_kind, -0.01300660_dbl_kind, & + 0.00988829_dbl_kind, -0.00336700_dbl_kind, 0.00333077_dbl_kind, & + -0.00235400_dbl_kind, 0.00140015_dbl_kind, 0.00100700_dbl_kind, & + 0.00085700_dbl_kind, 0.00064990_dbl_kind, 0.00059900_dbl_kind, & + 0.00037800_dbl_kind, -0.00033700_dbl_kind, 0.00027600_dbl_kind, & + 0.00018200_dbl_kind, -0.00017400_dbl_kind, -0.00012400_dbl_kind, & + 0.00001250_dbl_kind/) real (dbl_kind), parameter :: ecrate(pecclen) = & ! rates for eccen/fvelp cos/sin series - & (/ 4.2072050_dbl_kind, 7.3460910_dbl_kind, 17.8572630_dbl_kind, & - & 17.2205460_dbl_kind, 16.8467330_dbl_kind, 5.1990790_dbl_kind, & - & 18.2310760_dbl_kind, 26.2167580_dbl_kind, 6.3591690_dbl_kind, & - & 16.2100160_dbl_kind, 3.0651810_dbl_kind, 16.5838290_dbl_kind, & - & 18.4939800_dbl_kind, 6.1909530_dbl_kind, 18.8677930_dbl_kind, & - & 17.4255670_dbl_kind, 6.1860010_dbl_kind, 18.4174410_dbl_kind, & - & 0.6678630_dbl_kind/) + (/ 4.2072050_dbl_kind, 7.3460910_dbl_kind, 17.8572630_dbl_kind, & + 17.2205460_dbl_kind, 16.8467330_dbl_kind, 5.1990790_dbl_kind, & + 18.2310760_dbl_kind, 26.2167580_dbl_kind, 6.3591690_dbl_kind, & + 16.2100160_dbl_kind, 3.0651810_dbl_kind, 16.5838290_dbl_kind, & + 18.4939800_dbl_kind, 6.1909530_dbl_kind, 18.8677930_dbl_kind, & + 17.4255670_dbl_kind, 6.1860010_dbl_kind, 18.4174410_dbl_kind, & + 0.6678630_dbl_kind/) real (dbl_kind), parameter :: ecphas(pecclen) = & ! phases for eccen/fvelp cos/sin series - & (/ 28.620089_dbl_kind, 193.788772_dbl_kind, 308.307024_dbl_kind, & - & 320.199637_dbl_kind, 279.376984_dbl_kind, 87.195000_dbl_kind, & - & 349.129677_dbl_kind, 128.443387_dbl_kind, 154.143880_dbl_kind, & - & 291.269597_dbl_kind, 114.860583_dbl_kind, 332.092251_dbl_kind, & - & 296.414411_dbl_kind, 145.769910_dbl_kind, 337.237063_dbl_kind, & - & 152.092288_dbl_kind, 126.839891_dbl_kind, 210.667199_dbl_kind, & - & 72.108838_dbl_kind/) + (/ 28.620089_dbl_kind, 193.788772_dbl_kind, 308.307024_dbl_kind, & + 320.199637_dbl_kind, 279.376984_dbl_kind, 87.195000_dbl_kind, & + 349.129677_dbl_kind, 128.443387_dbl_kind, 154.143880_dbl_kind, & + 291.269597_dbl_kind, 114.860583_dbl_kind, 332.092251_dbl_kind, & + 296.414411_dbl_kind, 145.769910_dbl_kind, 337.237063_dbl_kind, & + 152.092288_dbl_kind, 126.839891_dbl_kind, 210.667199_dbl_kind, & + 72.108838_dbl_kind/) ! Sine series data for computation of moving vernal equinox longitude of ! perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees). real (dbl_kind), parameter :: mvamp (pmvelen) = & ! amplitudes for mvelp sine series - & (/ 7391.0225890_dbl_kind, 2555.1526947_dbl_kind, 2022.7629188_dbl_kind, & - & -1973.6517951_dbl_kind, 1240.2321818_dbl_kind, 953.8679112_dbl_kind, & - & -931.7537108_dbl_kind, 872.3795383_dbl_kind, 606.3544732_dbl_kind, & - & -496.0274038_dbl_kind, 456.9608039_dbl_kind, 346.9462320_dbl_kind, & - & -305.8412902_dbl_kind, 249.6173246_dbl_kind, -199.1027200_dbl_kind, & - & 191.0560889_dbl_kind, -175.2936572_dbl_kind, 165.9068833_dbl_kind, & - & 161.1285917_dbl_kind, 139.7878093_dbl_kind, -133.5228399_dbl_kind, & - & 117.0673811_dbl_kind, 104.6907281_dbl_kind, 95.3227476_dbl_kind, & - & 86.7824524_dbl_kind, 86.0857729_dbl_kind, 70.5893698_dbl_kind, & - & -69.9719343_dbl_kind, -62.5817473_dbl_kind, 61.5450059_dbl_kind, & - & -57.9364011_dbl_kind, 57.1899832_dbl_kind, -57.0236109_dbl_kind, & - & -54.2119253_dbl_kind, 53.2834147_dbl_kind, 52.1223575_dbl_kind, & - & -49.0059908_dbl_kind, -48.3118757_dbl_kind, -45.4191685_dbl_kind, & - & -42.2357920_dbl_kind, -34.7971099_dbl_kind, 34.4623613_dbl_kind, & - & -33.8356643_dbl_kind, 33.6689362_dbl_kind, -31.2521586_dbl_kind, & - & -30.8798701_dbl_kind, 28.4640769_dbl_kind, -27.1960802_dbl_kind, & - & 27.0860736_dbl_kind, -26.3437456_dbl_kind, 24.7253740_dbl_kind, & - & 24.6732126_dbl_kind, 24.4272733_dbl_kind, 24.0127327_dbl_kind, & - & 21.7150294_dbl_kind, -21.5375347_dbl_kind, 18.1148363_dbl_kind, & - & -16.9603104_dbl_kind, -16.1765215_dbl_kind, 15.5567653_dbl_kind, & - & 15.4846529_dbl_kind, 15.2150632_dbl_kind, 14.5047426_dbl_kind, & - & -14.3873316_dbl_kind, 13.1351419_dbl_kind, 12.8776311_dbl_kind, & - & 11.9867234_dbl_kind, 11.9385578_dbl_kind, 11.7030822_dbl_kind, & - & 11.6018181_dbl_kind, -11.2617293_dbl_kind, -10.4664199_dbl_kind, & - & 10.4333970_dbl_kind, -10.2377466_dbl_kind, 10.1934446_dbl_kind, & - & -10.1280191_dbl_kind, 10.0289441_dbl_kind, -10.0034259_dbl_kind/) + (/ 7391.0225890_dbl_kind, 2555.1526947_dbl_kind, 2022.7629188_dbl_kind, & + -1973.6517951_dbl_kind, 1240.2321818_dbl_kind, 953.8679112_dbl_kind, & + -931.7537108_dbl_kind, 872.3795383_dbl_kind, 606.3544732_dbl_kind, & + -496.0274038_dbl_kind, 456.9608039_dbl_kind, 346.9462320_dbl_kind, & + -305.8412902_dbl_kind, 249.6173246_dbl_kind, -199.1027200_dbl_kind, & + 191.0560889_dbl_kind, -175.2936572_dbl_kind, 165.9068833_dbl_kind, & + 161.1285917_dbl_kind, 139.7878093_dbl_kind, -133.5228399_dbl_kind, & + 117.0673811_dbl_kind, 104.6907281_dbl_kind, 95.3227476_dbl_kind, & + 86.7824524_dbl_kind, 86.0857729_dbl_kind, 70.5893698_dbl_kind, & + -69.9719343_dbl_kind, -62.5817473_dbl_kind, 61.5450059_dbl_kind, & + -57.9364011_dbl_kind, 57.1899832_dbl_kind, -57.0236109_dbl_kind, & + -54.2119253_dbl_kind, 53.2834147_dbl_kind, 52.1223575_dbl_kind, & + -49.0059908_dbl_kind, -48.3118757_dbl_kind, -45.4191685_dbl_kind, & + -42.2357920_dbl_kind, -34.7971099_dbl_kind, 34.4623613_dbl_kind, & + -33.8356643_dbl_kind, 33.6689362_dbl_kind, -31.2521586_dbl_kind, & + -30.8798701_dbl_kind, 28.4640769_dbl_kind, -27.1960802_dbl_kind, & + 27.0860736_dbl_kind, -26.3437456_dbl_kind, 24.7253740_dbl_kind, & + 24.6732126_dbl_kind, 24.4272733_dbl_kind, 24.0127327_dbl_kind, & + 21.7150294_dbl_kind, -21.5375347_dbl_kind, 18.1148363_dbl_kind, & + -16.9603104_dbl_kind, -16.1765215_dbl_kind, 15.5567653_dbl_kind, & + 15.4846529_dbl_kind, 15.2150632_dbl_kind, 14.5047426_dbl_kind, & + -14.3873316_dbl_kind, 13.1351419_dbl_kind, 12.8776311_dbl_kind, & + 11.9867234_dbl_kind, 11.9385578_dbl_kind, 11.7030822_dbl_kind, & + 11.6018181_dbl_kind, -11.2617293_dbl_kind, -10.4664199_dbl_kind, & + 10.4333970_dbl_kind, -10.2377466_dbl_kind, 10.1934446_dbl_kind, & + -10.1280191_dbl_kind, 10.0289441_dbl_kind, -10.0034259_dbl_kind/) real (dbl_kind), parameter :: mvrate(pmvelen) = & ! rates for mvelp sine series - & (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, & - & 0.636717_dbl_kind, 31.983787_dbl_kind, 3.138886_dbl_kind, & - & 30.973257_dbl_kind, 44.828336_dbl_kind, 0.991874_dbl_kind, & - & 0.373813_dbl_kind, 43.668246_dbl_kind, 32.246691_dbl_kind, & - & 30.599444_dbl_kind, 2.147012_dbl_kind, 10.511172_dbl_kind, & - & 42.681324_dbl_kind, 13.650058_dbl_kind, 0.986922_dbl_kind, & - & 9.874455_dbl_kind, 13.013341_dbl_kind, 0.262904_dbl_kind, & - & 0.004952_dbl_kind, 1.142024_dbl_kind, 63.219948_dbl_kind, & - & 0.205021_dbl_kind, 2.151964_dbl_kind, 64.230478_dbl_kind, & - & 43.836462_dbl_kind, 47.439436_dbl_kind, 1.384343_dbl_kind, & - & 7.437771_dbl_kind, 18.829299_dbl_kind, 9.500642_dbl_kind, & - & 0.431696_dbl_kind, 1.160090_dbl_kind, 55.782177_dbl_kind, & - & 12.639528_dbl_kind, 1.155138_dbl_kind, 0.168216_dbl_kind, & - & 1.647247_dbl_kind, 10.884985_dbl_kind, 5.610937_dbl_kind, & - & 12.658184_dbl_kind, 1.010530_dbl_kind, 1.983748_dbl_kind, & - & 14.023871_dbl_kind, 0.560178_dbl_kind, 1.273434_dbl_kind, & - & 12.021467_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, & - & 76.438310_dbl_kind, 4.280910_dbl_kind, 13.218362_dbl_kind, & - & 17.818769_dbl_kind, 8.359495_dbl_kind, 56.792707_dbl_kind, & - & 8.448301_dbl_kind, 1.978796_dbl_kind, 8.863925_dbl_kind, & - & 0.186365_dbl_kind, 8.996212_dbl_kind, 6.771027_dbl_kind, & - & 45.815258_dbl_kind, 12.002811_dbl_kind, 75.278220_dbl_kind, & - & 65.241008_dbl_kind, 18.870667_dbl_kind, 22.009553_dbl_kind, & - & 64.604291_dbl_kind, 11.498094_dbl_kind, 0.578834_dbl_kind, & - & 9.237738_dbl_kind, 49.747842_dbl_kind, 2.147012_dbl_kind, & - & 1.196895_dbl_kind, 2.133898_dbl_kind, 0.173168_dbl_kind/) + (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, & + 0.636717_dbl_kind, 31.983787_dbl_kind, 3.138886_dbl_kind, & + 30.973257_dbl_kind, 44.828336_dbl_kind, 0.991874_dbl_kind, & + 0.373813_dbl_kind, 43.668246_dbl_kind, 32.246691_dbl_kind, & + 30.599444_dbl_kind, 2.147012_dbl_kind, 10.511172_dbl_kind, & + 42.681324_dbl_kind, 13.650058_dbl_kind, 0.986922_dbl_kind, & + 9.874455_dbl_kind, 13.013341_dbl_kind, 0.262904_dbl_kind, & + 0.004952_dbl_kind, 1.142024_dbl_kind, 63.219948_dbl_kind, & + 0.205021_dbl_kind, 2.151964_dbl_kind, 64.230478_dbl_kind, & + 43.836462_dbl_kind, 47.439436_dbl_kind, 1.384343_dbl_kind, & + 7.437771_dbl_kind, 18.829299_dbl_kind, 9.500642_dbl_kind, & + 0.431696_dbl_kind, 1.160090_dbl_kind, 55.782177_dbl_kind, & + 12.639528_dbl_kind, 1.155138_dbl_kind, 0.168216_dbl_kind, & + 1.647247_dbl_kind, 10.884985_dbl_kind, 5.610937_dbl_kind, & + 12.658184_dbl_kind, 1.010530_dbl_kind, 1.983748_dbl_kind, & + 14.023871_dbl_kind, 0.560178_dbl_kind, 1.273434_dbl_kind, & + 12.021467_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, & + 76.438310_dbl_kind, 4.280910_dbl_kind, 13.218362_dbl_kind, & + 17.818769_dbl_kind, 8.359495_dbl_kind, 56.792707_dbl_kind, & + 8.448301_dbl_kind, 1.978796_dbl_kind, 8.863925_dbl_kind, & + 0.186365_dbl_kind, 8.996212_dbl_kind, 6.771027_dbl_kind, & + 45.815258_dbl_kind, 12.002811_dbl_kind, 75.278220_dbl_kind, & + 65.241008_dbl_kind, 18.870667_dbl_kind, 22.009553_dbl_kind, & + 64.604291_dbl_kind, 11.498094_dbl_kind, 0.578834_dbl_kind, & + 9.237738_dbl_kind, 49.747842_dbl_kind, 2.147012_dbl_kind, & + 1.196895_dbl_kind, 2.133898_dbl_kind, 0.173168_dbl_kind/) real (dbl_kind), parameter :: mvphas(pmvelen) = & ! phases for mvelp sine series - & (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, & - & 348.1074_dbl_kind, 292.7252_dbl_kind, 165.1686_dbl_kind, & - & 263.7951_dbl_kind, 15.3747_dbl_kind, 58.5749_dbl_kind, & - & 40.8226_dbl_kind, 308.4258_dbl_kind, 240.0099_dbl_kind, & - & 222.9725_dbl_kind, 106.5937_dbl_kind, 114.5182_dbl_kind, & - & 268.7809_dbl_kind, 279.6869_dbl_kind, 39.6448_dbl_kind, & - & 126.4108_dbl_kind, 291.5795_dbl_kind, 307.2848_dbl_kind, & - & 18.9300_dbl_kind, 273.7596_dbl_kind, 143.8050_dbl_kind, & - & 191.8927_dbl_kind, 125.5237_dbl_kind, 172.7351_dbl_kind, & - & 316.7998_dbl_kind, 319.6024_dbl_kind, 69.7526_dbl_kind, & - & 123.5968_dbl_kind, 217.6432_dbl_kind, 85.5882_dbl_kind, & - & 156.2147_dbl_kind, 66.9489_dbl_kind, 20.2082_dbl_kind, & - & 250.7568_dbl_kind, 48.0188_dbl_kind, 8.3739_dbl_kind, & - & 17.0374_dbl_kind, 155.3409_dbl_kind, 94.1709_dbl_kind, & - & 221.1120_dbl_kind, 28.9300_dbl_kind, 117.1498_dbl_kind, & - & 320.5095_dbl_kind, 262.3602_dbl_kind, 336.2148_dbl_kind, & - & 233.0046_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, & - & 267.2772_dbl_kind, 78.9281_dbl_kind, 123.4722_dbl_kind, & - & 188.7132_dbl_kind, 180.1364_dbl_kind, 49.1382_dbl_kind, & - & 152.5268_dbl_kind, 98.2198_dbl_kind, 97.4808_dbl_kind, & - & 221.5376_dbl_kind, 168.2438_dbl_kind, 161.1199_dbl_kind, & - & 55.0196_dbl_kind, 262.6495_dbl_kind, 200.3284_dbl_kind, & - & 201.6651_dbl_kind, 294.6547_dbl_kind, 99.8233_dbl_kind, & - & 213.5577_dbl_kind, 154.1631_dbl_kind, 232.7153_dbl_kind, & - & 138.3034_dbl_kind, 204.6609_dbl_kind, 106.5938_dbl_kind, & - & 250.4676_dbl_kind, 332.3345_dbl_kind, 27.3039_dbl_kind/) + (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, & + 348.1074_dbl_kind, 292.7252_dbl_kind, 165.1686_dbl_kind, & + 263.7951_dbl_kind, 15.3747_dbl_kind, 58.5749_dbl_kind, & + 40.8226_dbl_kind, 308.4258_dbl_kind, 240.0099_dbl_kind, & + 222.9725_dbl_kind, 106.5937_dbl_kind, 114.5182_dbl_kind, & + 268.7809_dbl_kind, 279.6869_dbl_kind, 39.6448_dbl_kind, & + 126.4108_dbl_kind, 291.5795_dbl_kind, 307.2848_dbl_kind, & + 18.9300_dbl_kind, 273.7596_dbl_kind, 143.8050_dbl_kind, & + 191.8927_dbl_kind, 125.5237_dbl_kind, 172.7351_dbl_kind, & + 316.7998_dbl_kind, 319.6024_dbl_kind, 69.7526_dbl_kind, & + 123.5968_dbl_kind, 217.6432_dbl_kind, 85.5882_dbl_kind, & + 156.2147_dbl_kind, 66.9489_dbl_kind, 20.2082_dbl_kind, & + 250.7568_dbl_kind, 48.0188_dbl_kind, 8.3739_dbl_kind, & + 17.0374_dbl_kind, 155.3409_dbl_kind, 94.1709_dbl_kind, & + 221.1120_dbl_kind, 28.9300_dbl_kind, 117.1498_dbl_kind, & + 320.5095_dbl_kind, 262.3602_dbl_kind, 336.2148_dbl_kind, & + 233.0046_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, & + 267.2772_dbl_kind, 78.9281_dbl_kind, 123.4722_dbl_kind, & + 188.7132_dbl_kind, 180.1364_dbl_kind, 49.1382_dbl_kind, & + 152.5268_dbl_kind, 98.2198_dbl_kind, 97.4808_dbl_kind, & + 221.5376_dbl_kind, 168.2438_dbl_kind, 161.1199_dbl_kind, & + 55.0196_dbl_kind, 262.6495_dbl_kind, 200.3284_dbl_kind, & + 201.6651_dbl_kind, 294.6547_dbl_kind, 99.8233_dbl_kind, & + 213.5577_dbl_kind, 154.1631_dbl_kind, 232.7153_dbl_kind, & + 138.3034_dbl_kind, 204.6609_dbl_kind, 106.5938_dbl_kind, & + 250.4676_dbl_kind, 332.3345_dbl_kind, 27.3039_dbl_kind/) !---------------------------Local variables---------------------------------- integer(int_kind) :: i ! Index for series summations @@ -571,7 +571,7 @@ SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & obsum = 0.0_dbl_kind do i = 1, poblen obsum = obsum + obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + & - & obphas(i))*degrad) + obphas(i))*degrad) end do obliq = 23.320556_dbl_kind + obsum @@ -629,7 +629,7 @@ SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & mvsum = 0.0_dbl_kind do i = 1, pmvelen mvsum = mvsum + mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + & - & mvphas(i))*degrad) + mvphas(i))*degrad) end do mvelp = fvelp/degrad + 50.439273_dbl_kind*psecdeg*years + 3.392506_dbl_kind + mvsum @@ -670,8 +670,8 @@ SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & ! 1978) is 0. lambm0 = 2._dbl_kind*((.5_dbl_kind*eccen + .125_dbl_kind*eccen3)*(1._dbl_kind + beta)*sin(mvelpp) & - & - .250_dbl_kind*eccen2*(.5_dbl_kind + beta)*sin(2._dbl_kind*mvelpp) & - & + .125_dbl_kind*eccen3*(1._dbl_kind/3._dbl_kind + beta)*sin(3._dbl_kind*mvelpp)) + - .250_dbl_kind*eccen2*(.5_dbl_kind + beta)*sin(2._dbl_kind*mvelpp) & + + .125_dbl_kind*eccen3*(1._dbl_kind/3._dbl_kind + beta)*sin(3._dbl_kind*mvelpp)) if ( log_print ) then write(warnstr,F03) subname//'------ Computed Orbital Parameters ------' @@ -759,7 +759,7 @@ SUBROUTINE shr_orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr ,delta ,eccf) sinl = sin(lmm) lamb = lambm + eccen*(2._dbl_kind*sinl + eccen*(1.25_dbl_kind*sin(2._dbl_kind*lmm) & - & + eccen*((13.0_dbl_kind/12.0_dbl_kind)*sin(3._dbl_kind*lmm) - 0.25_dbl_kind*sinl))) + + eccen*((13.0_dbl_kind/12.0_dbl_kind)*sin(3._dbl_kind*lmm) - 0.25_dbl_kind*sinl))) ! Using the obliquity, eccentricity, moving vernal equinox longitude of ! perihelion (plus), and earths true longitude, the declination (delta) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 2240f40a9..9033e487b 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -4323,6 +4323,8 @@ subroutine icepack_step_radiation (dt, ncat, & end subroutine icepack_step_radiation +!======================================================================= + ! Delta-Eddington solution expressions !======================================================================= diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index ba2bc8134..c0ddc965e 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -305,11 +305,24 @@ subroutine linear_itd (ncat, hin_max, & if (hicen_init(n) > puny .and. & hicen_init(n+1) > puny) then - ! interpolate between adjacent category growth rates - slope = (dhicen(n+1) - dhicen(n)) / & + + if ((hicen_init(n+1) - hicen_init(n))>0) then + + ! interpolate between adjacent category growth rates + slope = (dhicen(n+1) - dhicen(n)) / & (hicen_init(n+1) - hicen_init(n)) - hbnew(n) = hin_max(n) + dhicen(n) & + hbnew(n) = hin_max(n) + dhicen(n) & + slope * (hin_max(n) - hicen_init(n)) + + else + + write(warnstr,*) subname, & + 'ITD Thermodynamics: hicen_init(n+1) <= hicen_init(n)' + call icepack_warnings_add(warnstr) + call icepack_warnings_setabort(.false.) + + endif + elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0 hbnew(n) = hin_max(n) + dhicen(n) elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0 diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index 1a28ed1c5..c50aff2db 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -217,32 +217,28 @@ subroutine temperature_changes_salinity(dt, & !----------------------------------------------------------------- !mclaren: Should there be an if calc_Tsfc statement here then?? + dswabs = c0 if (sw_redist) then + dt_rhoi_hlyr = dt / (rhoi*hilyr) + do k = 1, nilyr + Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc + Tmlt = liquidus_temperature_mush(zSin(k)) - dt_rhoi_hlyr = dt / (rhoi*hilyr) - - do k = 1, nilyr - - Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc - - Tmlt = liquidus_temperature_mush(zSin(k)) - - if (zTin(k) <= Tmlt - sw_dtemp) then - ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2) - Iswabs_tmp = min(Iswabs(k), & - sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr) - endif - if (Iswabs_tmp < puny) Iswabs_tmp = c0 - - dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) - - fswsfc = fswsfc + dswabs - fswint = fswint - dswabs - Iswabs(k) = Iswabs_tmp - - enddo - + if (zTin(k) <= Tmlt - sw_dtemp) then + ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2) + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr) + endif + if (Iswabs_tmp < puny) Iswabs_tmp = c0 + dswabs = dswabs + min(Iswabs(k) - Iswabs_tmp, fswint) + Iswabs(k) = Iswabs_tmp + enddo + endif + if (.not. lsnow) then ! hs <= hs_min + dswabs = dswabs + sum(Sswabs(:)) endif + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs if (lsnow) then ! case with snow diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 74cc7cdc2..45db195d9 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1887,7 +1887,7 @@ subroutine conservation_check_vthermo(dt, & - fsnow*Lfresh - fadvocn) * dt ferr = abs(efinal-einit-einp) / dt - if (ferr > ferrmax) then + if (ferr > 1.1_dbl_kind*ferrmax) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" ) @@ -2111,7 +2111,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & dsnown , & lmask_n , lmask_s , & mlt_onset , frz_onset , & - yday , prescribed_ice) + yday , prescribed_ice, & + zlvs) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -2137,7 +2138,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & aice , & ! sea ice concentration vice , & ! volume per unit area of ice (m) vsno , & ! volume per unit area of snow (m) - zlvl , & ! atm level height (m) + zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) uatm , & ! wind velocity components (m/s) vatm , & wind , & ! wind speed (m/s) @@ -2221,7 +2222,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & real (kind=dbl_kind), optional, intent(in) :: & HDO_ocn , & ! ocean concentration of HDO (kg/kg) H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) - H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) + H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) + zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice @@ -2539,7 +2541,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & Qa_iso=l_Qa_iso, & Qref_iso=Qrefn_iso, & uvel=uvel, vvel=vvel, & - Uref=Urefn) + Uref=Urefn, zlvs=zlvs) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc or calc_strair diff --git a/columnphysics/version.txt b/columnphysics/version.txt index 1792f6964..9d36f27cb 100644 --- a/columnphysics/version.txt +++ b/columnphysics/version.txt @@ -1 +1 @@ -ICEPACK 1.2.4 +ICEPACK 1.2.5 diff --git a/configuration/scripts/icepack.batch.csh b/configuration/scripts/icepack.batch.csh index 1549ba197..52b280ed3 100755 --- a/configuration/scripts/icepack.batch.csh +++ b/configuration/scripts/icepack.batch.csh @@ -112,6 +112,21 @@ cat >> ${jobfile} << EOFB ###SBATCH --mail-user username@domain.com EOFB +else if (${ICE_MACHINE} =~ compy*) then +cat >> ${jobfile} << EOFB +#SBATCH -J ${ICE_CASENAME} +#SBATCH -A ${acct} +#SBATCH --qos ${ICE_MACHINE_QUEUE} +#SBATCH --time ${ICE_RUNLENGTH} +#SBATCH --nodes ${nnodes} +#SBATCH --ntasks ${ntasks} +#SBATCH --cpus-per-task ${nthrds} +###SBATCH -e filename +###SBATCH -o filename +###SBATCH --mail-type FAIL +###SBATCH --mail-user username@domain.com +EOFB + else if (${ICE_MACHINE} =~ badger*) then cat >> ${jobfile} << EOFB #SBATCH -J ${ICE_CASENAME} @@ -125,6 +140,15 @@ cat >> ${jobfile} << EOFB #SBATCH --qos=standby EOFB +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=${ICE_RUNLENGTH} +#PBS -W umask=022 +EOFB + else if (${ICE_MACHINE} =~ high_Sierra*) then cat >> ${jobfile} << EOFB # nothing to do diff --git a/configuration/scripts/icepack.launch.csh b/configuration/scripts/icepack.launch.csh index 7acd7054b..dc7ddd489 100755 --- a/configuration/scripts/icepack.launch.csh +++ b/configuration/scripts/icepack.launch.csh @@ -9,7 +9,7 @@ set jobfile = $1 #========================================== -if (${ICE_MACHINE} =~ gordon* || ${ICE_MACHINE} =~ conrad* || ${ICE_MACHINE} =~ onyx* ) then +if (${ICE_MACHINE} =~ gordon* || ${ICE_MACHINE} =~ conrad* || ${ICE_MACHINE} =~ onyx* || ${ICE_MACHINE} =~ daley* || ${ICE_MACHINE} =~ banting* ) then cat >> ${jobfile} << EOFR aprun -n 1 -N 1 -d 1 ./icepack >&! \$ICE_RUNLOG_FILE EOFR diff --git a/configuration/scripts/machines/Macros.banting_gnu b/configuration/scripts/machines/Macros.banting_gnu new file mode 100644 index 000000000..3652712b1 --- /dev/null +++ b/configuration/scripts/machines/Macros.banting_gnu @@ -0,0 +1,31 @@ +#============================================================================== +# Makefile macros for ECCC banting +#============================================================================== +# For use with GNU compiler +#============================================================================== + +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 +endif + +SCC := cc +SFC := ftn +FC := $(SFC) +CC := $(SCC) +LD:= $(FC) + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -fopenmp + CFLAGS += -fopenmp + FFLAGS += -fopenmp +endif diff --git a/configuration/scripts/machines/Macros.banting_intel b/configuration/scripts/machines/Macros.banting_intel new file mode 100644 index 000000000..a87a58769 --- /dev/null +++ b/configuration/scripts/machines/Macros.banting_intel @@ -0,0 +1,33 @@ +#============================================================================== +# Makefile macros for ECCC banting +#============================================================================== +# For use with intel compiler +#============================================================================== + +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise +#-xHost + +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 +#-init=snan,arrays +# -heap-arrays 1024 +else + FFLAGS += -O2 +endif + +SCC := cc +SFC := ftn +FC := $(SFC) +CC := $(SCC) +LD:= $(FC) + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp +endif diff --git a/configuration/scripts/machines/Macros.compy_intel b/configuration/scripts/machines/Macros.compy_intel new file mode 100644 index 000000000..a65e93975 --- /dev/null +++ b/configuration/scripts/machines/Macros.compy_intel @@ -0,0 +1,39 @@ +#============================================================================== +# Makefile macro for PNNL compy, intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise -xHost + +FIXEDFLAGS := -132 +FREEFLAGS := -FR +FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -xHost +FFLAGS_NOOPT:= -O0 + +ifeq ($(ICE_BLDDEBUG), true) + FFLAGS += -O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created +else + FFLAGS += -O2 +endif + +SCC := icc +SFC := ifort +CC := $(SCC) +FC := $(SFC) +LD := $(FC) + +#NETCDF_PATH := /share/apps/netcdf/4.6.3/intel/19.0.3 +INCLDIR := $(INCLDIR) + +#LIB_NETCDF := $(NETCDF_PATH)/lib +#LIB_PNETCDF := $(PNETCDF_PATH)/lib +#LIB_MPI := $(IMPILIBDIR) +#SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -openmp + CFLAGS += -openmp + FFLAGS += -openmp +endif + diff --git a/configuration/scripts/machines/Macros.conda_macos b/configuration/scripts/machines/Macros.conda_macos index 6f5274b61..583a37773 100644 --- a/configuration/scripts/machines/Macros.conda_macos +++ b/configuration/scripts/machines/Macros.conda_macos @@ -32,6 +32,7 @@ ifeq ($(strip $(SDKPATH)),) CFLAGS_HOST := else CFLAGS_HOST = -isysroot $(SDKPATH) + LD += -L$(SDKPATH)/usr/lib endif # Necessary flag to compile with OpenMP support diff --git a/configuration/scripts/machines/Macros.daley_gnu b/configuration/scripts/machines/Macros.daley_gnu new file mode 100644 index 000000000..8daf17fba --- /dev/null +++ b/configuration/scripts/machines/Macros.daley_gnu @@ -0,0 +1,31 @@ +#============================================================================== +# Makefile macros for ECCC daley +#============================================================================== +# For use with GNU compiler +#============================================================================== + +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 +endif + +SCC := cc +SFC := ftn +FC := $(SFC) +CC := $(SCC) +LD:= $(FC) + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -fopenmp + CFLAGS += -fopenmp + FFLAGS += -fopenmp +endif diff --git a/configuration/scripts/machines/Macros.daley_intel b/configuration/scripts/machines/Macros.daley_intel new file mode 100644 index 000000000..2c0515b05 --- /dev/null +++ b/configuration/scripts/machines/Macros.daley_intel @@ -0,0 +1,33 @@ +#============================================================================== +# Makefile macros for ECCC daley +#============================================================================== +# For use with intel compiler +#============================================================================== + +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise +#-xHost + +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 +# -init=snan,arrays +# -heap-arrays 1024 +else + FFLAGS += -O2 +endif + +SCC := cc +SFC := ftn +FC := $(SFC) +CC := $(SCC) +LD:= $(FC) + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp +endif diff --git a/configuration/scripts/machines/env.banting_gnu b/configuration/scripts/machines/env.banting_gnu new file mode 100755 index 000000000..4ce2be5bb --- /dev/null +++ b/configuration/scripts/machines/env.banting_gnu @@ -0,0 +1,35 @@ +#!/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 unload craype-x86-skylake +module unload PrgEnv-intel + +module load PrgEnv-gnu # GNU compiler +module load craype-x86-skylake # Reload +module load cray-netcdf # NetCDF +module load cray-hdf5 # HDF5 +setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME banting +setenv ICE_MACHINE_ENVNAME gnu +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR ~/data/banting/icepack/runs +setenv ICE_MACHINE_INPUTDATA /home/ords/cmdd/cmde/sice500/ +setenv ICE_MACHINE_BASELINE ~/data/banting/icepack/baselines +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 40 +setenv ICE_MACHINE_MAXRUNLENGTH 3 +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.banting_intel b/configuration/scripts/machines/env.banting_intel new file mode 100755 index 000000000..08da69cf9 --- /dev/null +++ b/configuration/scripts/machines/env.banting_intel @@ -0,0 +1,32 @@ +#!/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 +module load PrgEnv-intel # Intel compiler +module load cray-netcdf # NetCDF +module load cray-hdf5 # HDF5 +setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME banting +setenv ICE_MACHINE_MACHINFO "Cray XC50, Intel Xeon Gold 6148 (Skylake)" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO "Intel 19.0.3.199, cray-mpich/7.7.7, cray-netcdf/4.6.1.3" +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR ~/data/banting/icepack/runs +setenv ICE_MACHINE_INPUTDATA /home/ords/cmdd/cmde/sice500/ +setenv ICE_MACHINE_BASELINE ~/data/banting/icepack/baselines +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 40 +setenv ICE_MACHINE_MAXRUNLENGTH 3 +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.compy_intel b/configuration/scripts/machines/env.compy_intel new file mode 100755 index 000000000..378249423 --- /dev/null +++ b/configuration/scripts/machines/env.compy_intel @@ -0,0 +1,39 @@ +#!/bin/csh -f + +set inp = "undefined" +if ($#argv == 1) then + set inp = $1 +endif + +if ("$inp" != "-nomodules") then + +source /share/apps/modules/init/csh + +module purge +module load intel/19.0.5 +#module load intelmpi/2019u4 +module load netcdf/4.6.3 +module load hdf5/1.10.5 + +#setenv NETCDF_PATH ${NETCDF_DIR} +setenv NETCDF_PATH /share/apps/netcdf/4.6.3/intel/19.0.5 +limit coredumpsize unlimited +limit stacksize unlimited + +endif + +setenv ICE_MACHINE_MACHNAME compy +setenv ICE_MACHINE_MACHINFO "PNNL Intel Xeon Skylake with 192 GB of DDR4 DRAM" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO "intel/19.0.5" +setenv ICE_MACHINE_MAKE gmake +setenv ICE_MACHINE_WKDIR /compyfs/$USER/ICEPACK_RUNS +setenv ICE_MACHINE_INPUTDATA /compyfs/inputdata/cice-consortium/ +setenv ICE_MACHINE_BASELINE /compyfs/$USER/ICEPACK_BASELINE +setenv ICE_MACHINE_SUBMIT "sbatch " +setenv ICE_MACHINE_ACCT e3sm +setenv ICE_MACHINE_QUEUE "short" +setenv ICE_MACHINE_TPNODE 40 # tasks per node +setenv ICE_MACHINE_BLDTHRDS 4 +setenv ICE_MACHINE_QSTAT "squeue --jobs=" + diff --git a/configuration/scripts/machines/env.daley_gnu b/configuration/scripts/machines/env.daley_gnu new file mode 100755 index 000000000..a091ea706 --- /dev/null +++ b/configuration/scripts/machines/env.daley_gnu @@ -0,0 +1,35 @@ +#!/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 unload craype-x86-skylake +module unload PrgEnv-intel + +module load PrgEnv-gnu # GNU compiler +module load craype-x86-skylake # Reload +module load cray-netcdf # NetCDF +module load cray-hdf5 # HDF5 +setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME daley +setenv ICE_MACHINE_ENVNAME gnu +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR ~/data/daley/icepack/runs +setenv ICE_MACHINE_INPUTDATA /home/ords/cmdd/cmde/sice500/ +setenv ICE_MACHINE_BASELINE ~/data/daley/icepack/baselines +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 40 +setenv ICE_MACHINE_MAXRUNLENGTH 3 +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.daley_intel b/configuration/scripts/machines/env.daley_intel new file mode 100755 index 000000000..3f0a3bf7d --- /dev/null +++ b/configuration/scripts/machines/env.daley_intel @@ -0,0 +1,32 @@ +#!/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 +module load PrgEnv-intel # Intel compiler +module load cray-netcdf # NetCDF +module load cray-hdf5 # HDF5 +setenv HDF5_USE_FILE_LOCKING FALSE # necessary since data is on an NFS filesystem + +endif + +setenv ICE_MACHINE_MACHNAME daley +setenv ICE_MACHINE_MACHINFO "Cray XC50, Intel Xeon Gold 6148 (Skylake)" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO "Intel 19.0.3.199, cray-mpich/7.7.6, cray-netcdf/4.6.1.3" +setenv ICE_MACHINE_MAKE make +setenv ICE_MACHINE_WKDIR ~/data/daley/icepack/runs +setenv ICE_MACHINE_INPUTDATA /home/ords/cmdd/cmde/sice500/ +setenv ICE_MACHINE_BASELINE ~/data/daley/icepack/baselines +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_TPNODE 40 +setenv ICE_MACHINE_MAXRUNLENGTH 3 +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/tests/report_results.csh b/configuration/scripts/tests/report_results.csh index 69ab702be..df53cd5e2 100755 --- a/configuration/scripts/tests/report_results.csh +++ b/configuration/scripts/tests/report_results.csh @@ -82,11 +82,11 @@ if ("${shrepo}" !~ "*cice-consortium*") then endif set noglob -set green = "\![#00C000](https://placehold.it/15/00C000/000000?text=+)" -set red = "\![#F00000](https://placehold.it/15/F00000/000000?text=+)" -set orange = "\![#FFA500](https://placehold.it/15/FFA500/000000?text=+)" -set yellow = "\![#FFE600](https://placehold.it/15/FFE600/000000?text=+)" -set gray = "\![#AAAAAA](https://placehold.it/15/AAAAAA/000000?text=+)" +set green = "\![#00C000](images/00C000.png)" +set red = "\![#F00000](images/F00000.png)" +set orange = "\![#FFA500](images/FFA500.png)" +set yellow = "\![#FFE600](images/FFE600.png)" +set gray = "\![#AAAAAA](images/AAAAAA.png)" unset noglob #============================================================== diff --git a/doc/source/conf.py b/doc/source/conf.py index 86e179556..37dc338bd 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -62,9 +62,9 @@ # built documents. # # The short X.Y version. -version = u'1.2.4' +version = u'1.2.5' # The full version, including alpha/beta/rc tags. -version = u'1.2.4' +version = u'1.2.5' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/developer_guide/dg_adding_diagnostics.rst b/doc/source/developer_guide/dg_adding_diagnostics.rst new file mode 100755 index 000000000..ef5960830 --- /dev/null +++ b/doc/source/developer_guide/dg_adding_diagnostics.rst @@ -0,0 +1,63 @@ +:tocdepth: 3 + +.. _adddiag: + +Adding diagnostics +================== + +Icepack produces separate ASCII (text) log output for four cells, each with a different initial condition (full ITD, slab ice, ice free, land) designated by the variable ``n`` here. Each of the diagnostic files contains the state information for that cell. The procedure for adding diagnostic variables to the output is outlined here. + +#. For non-BGC variables, edit **icedrv\_diagnostics.F90**: + + - If the variable is already defined within the code, then add it to a "use" statement in the subroutine + ``runtime_diags``. + + - Note that if the variable is not readily accessible through a use statement, then a global variable may need to + be defined. This might be in **icedrv\_state.F90** or **icedrv\_flux.F90** for example. + + - Additionally, if the variable is a derived quantity, then the variables needed to calculate the new quantity + may need to be added to a use statement. For example, see how ``hiavg`` and ``hsavg`` are computed. + + - If the variable is a scalar, then follow the example of ``aice`` or ``hiavg``, copying the write statement to + an appropriate place in the output list, and editing as needed. The format "900" is appropriate for most scalars. + The following example adds snow melt (``melts``). + + .. code-block:: fortran + + use icedrv_flux, only: melts + + write(nu_diag_out+n-1,900) 'snow melt (m) = ',melts(n) ! snow melt + + - If the variable is an array, then you can compute the mean value (e.g. ``hiavg``) or print the array values (e.g. ``fiso_evap``). + This may requires adding the array sizes and a counter for the loop(s). E.g. to print + the category ice area, ``aicen`` over ncat thickness categories: + + .. code-block:: fortran + + use icedrv_domain_size, only: ncat + + use icedrv_state, only: aicen + + ! local variables + + integer (kind=int_kind) :: & + n, nc, k + + do nc = 1,ncat + write(nu_diag_out+n-1,901) 'Category ice area = ',aicen(n,nc),nc ! category ice area + enddo + + - If the variable is a tracer, then in addition to the variable trcr or trcrn, you will need the tracer + index (e.g. ``nt_Tsfc``). + + - In some cases, a new format statement might be needed. + +#. For BGC variables, edit **icedrv\_diagnostics\_bgc.F90**: + + - If the variable is already defined within the code, then add it to a "use" statement in the subroutine + ``hbrine_diags`` or ``bgc_diags`` or ``zsal_diags`` and follow a similar procedure for state variables as above. + + - Note that the BGC needs to be activated and the particular tracer turned on. + +In general, try to format the output statements to line up with the surrounding print messages. This may require a couple of tries to get it to compile and run. + diff --git a/doc/source/developer_guide/index.rst b/doc/source/developer_guide/index.rst index 31757e2a6..8206dbd3d 100755 --- a/doc/source/developer_guide/index.rst +++ b/doc/source/developer_guide/index.rst @@ -16,4 +16,5 @@ Developer Guide dg_driver.rst dg_scripts.rst dg_adding_tracers.rst + dg_adding_diagnostics.rst dg_other.rst diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index cc3a98356..e1c1816fd 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -509,7 +509,8 @@ either Celsius or Kelvin units). "yday", "day of the year", "" "year_init", ":math:`\bullet` the initial year", "" "**Z**", "", "" - "zlvl", "atmospheric level height", "m" + "zlvl", "atmospheric level height for momentum (and scalars if zlvs not present)", "m" + "zlvs", "atmospheric level height for scalars", "m" "zref", "reference height for stability", "10. m" "zTrf", "reference height for :math:`T_{ref}`, :math:`Q_{ref}`, :math:`U_{ref}`", "2. m" "zvir", "gas constant (water vapor)/gas constant (air) - 1", "0.606" diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index a763b1fa5..1a1e1b73e 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -579,8 +579,8 @@ @Article{Jeffery14 url = {http://dx.doi.org/10.1002/2013JC009634} } @Article{Saha14 - author = "S. Saha and S. Moorthi and X. Wu and J. Wang and Coauthors", - title = "{Two modes of sea-ice gravity drainage: a parameterization for large-scale modeling}", + author = "S. Saha and S. Moorthi and X. Wu and J. Wang and S. Nadiga and P. Tripp and D Behringer and Y. Hou and H. Chuang and M. Iredell and M. Ek and J. Meng and R. Yang and M.P. Mendez and H. van den Dool and Q. Zhang and W. Wang and M. Chen and E. Becker", + title = "{The NCEP Climate Forecast System Version 2}", journal = JC, year = {2014}, volume = {27}, diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst index 2e3507f90..763a757d0 100755 --- a/doc/source/science_guide/sg_boundary_forcing.rst +++ b/doc/source/science_guide/sg_boundary_forcing.rst @@ -12,6 +12,7 @@ Atmosphere and ocean boundary forcing :widths: 10, 25, 25 ":math:`z_o`", "Atmosphere level height", "From *atmosphere model* to *sea ice model*" + ":math:`z_{o,s}`", "Atmosphere level height (scalar quantities) (optional)", "From *atmosphere model* to *sea ice model*" ":math:`\vec{U}_a`", "Wind velocity", "From *atmosphere model* to *sea ice model*" ":math:`Q_a`", "Specific humidity", "From *atmosphere model* to *sea ice model*" ":math:`\rho_a`", "Air density", "From *atmosphere model* to *sea ice model*" @@ -94,7 +95,7 @@ Atmosphere ---------- The wind velocity, specific humidity, air density and potential -temperature at the given level height :math:`z_\circ` are used to +temperature at the given level height :math:`z_\circ` (optionally :math:`z_{\circ,s}`, see below) are used to compute transfer coefficients used in formulas for the surface wind stress and turbulent heat fluxes :math:`\vec\tau_a`, :math:`F_s`, and :math:`F_l`, as described below. The sensible and latent heat fluxes, @@ -121,7 +122,7 @@ turbulent heat fluxes are computed in subroutine The resulting equations are provided here. The wind stress and turbulent heat flux calculation accounts for both -stable and unstable atmosphere–ice boundary layers. Define the +stable and unstable atmosphere–ice boundary layers. We first define the "stability" .. math:: @@ -133,7 +134,7 @@ stable and unstable atmosphere–ice boundary layers. Define the where :math:`\kappa` is the von Karman constant, :math:`g` is gravitational acceleration, and :math:`u^*`, :math:`\Theta^*` and :math:`Q^*` are turbulent scales for velocity difference, temperature, and humidity, -respectively, given the ice velocity :math:`\vec{U}_i`: +respectively, computed as (given the ice velocity :math:`\vec{U}_i`): .. math:: \begin{aligned} @@ -143,6 +144,10 @@ respectively, given the ice velocity :math:`\vec{U}_i`: \end{aligned} :label: stars + +Note that *atmo_boundary_layer* also accepts an optional argument, ``zlvs``, to support staggered atmospheric levels, i.e. receiving scalar quantities from the atmospheric model (humidity and temperature) +at a different vertical level than the winds. In that case a separate stability :math:`\Upsilon_s` is computed using the same formula as above but substituting :math:`z_o` by :math:`z_{o,s}`. + Within the :math:`u^*` expression, :math:`U_{\Delta\textrm{min}}` is the minimum allowable value of :math:`|\vec{U}_{a} - \vec{U}_{i}|` , which is set to of 0.5 m/s for high frequency coupling (``highfreq`` =.true.). When high frequency coupling is turned off (``highfreq`` =.false.), it is assumed in equation :eq:`stars` that: @@ -197,11 +202,11 @@ The coefficients are then updated as .. math:: \begin{aligned} c_u^\prime&=&{c_u\over 1+c_u\left(\lambda-\psi_m\right)/\kappa} \\ - c_\theta^\prime&=& {c_\theta\over 1+c_\theta\left(\lambda-\psi_s\right)/\kappa}\\ + c_\theta^\prime&=& {c_\theta\over 1+c_\theta\left(\lambda_s-\psi_s\right)/\kappa}\\ c_q^\prime&=&c_\theta^\prime\end{aligned} :label: coeff1 -where :math:`\lambda = \ln\left(z_\circ/z_{ref}\right)`. The +where :math:`\lambda = \ln\left(z_\circ/z_{ref}\right)` and :math:`\lambda_s = \ln\left(z_{\circ,s}/z_{ref}\right)` if staggered atmospheric levels are used, else :math:`\lambda_s=\lambda`. The first iteration ends with new turbulent scales from equations :eq:`stars`. After ``natmiter`` iterations the latent and sensible heat flux coefficients are computed, along with the wind stress: diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index ad942f320..aff4f9ccc 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -743,9 +743,9 @@ Climate System Model (CCSM3), the albedo depends on the temperature and thickness of ice and snow and on the spectral distribution of the incoming solar radiation. Albedo parameters have been chosen to fit observations from the SHEBA field experiment. For -:math:`T_{sf} < -1^{\circ}C` and :math:`h_i > `\ ``ahmax``, the bare ice -albedo is 0.78 for visible wavelengths (:math:`<700`\ nm) and 0.36 for -near IR wavelengths (:math:`>700`\ nm). As :math:`h_i` decreases from +:math:`T_{sf} < -1^{\circ}C` and :math:`h_i >` \ `ahmax`, the bare ice +albedo is 0.78 for visible wavelengths (:math:`<700` \ nm) and 0.36 for +near IR wavelengths (:math:`>700` \ nm). As :math:`h_i` decreases from ahmax to zero, the ice albedo decreases smoothly (using an arctangent function) to the ocean albedo, 0.06. The ice albedo in both spectral bands decreases by 0.075 as :math:`T_{sf}` rises from diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include index e55e989cf..bb9f65556 100644 --- a/doc/source/user_guide/interfaces.include +++ b/doc/source/user_guide/interfaces.include @@ -23,7 +23,7 @@ icepack_atm_boundary Cdn_atm_ratio_n, & Qa_iso, Qref_iso, & uvel, vvel, & - Uref) + Uref, zlvs) character (len=3), intent(in) :: & sfctype ! ice or ocean @@ -34,7 +34,7 @@ icepack_atm_boundary uatm , & ! x-direction wind speed (m/s) vatm , & ! y-direction wind speed (m/s) wind , & ! wind speed (m/s) - zlvl , & ! atm level height (m) + zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) @@ -63,7 +63,8 @@ icepack_atm_boundary real (kind=dbl_kind), optional, intent(in) :: & uvel , & ! x-direction ice speed (m/s) - vvel ! y-direction ice speed (m/s) + vvel , & ! y-direction ice speed (m/s) + zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), optional, intent(out) :: & Uref ! reference height wind speed (m/s) @@ -2067,7 +2068,8 @@ icepack_step_therm1 dsnown , & lmask_n , lmask_s , & mlt_onset , frz_onset , & - yday , prescribed_ice) + yday , prescribed_ice, & + zlvs) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -2093,7 +2095,7 @@ icepack_step_therm1 aice , & ! sea ice concentration vice , & ! volume per unit area of ice (m) vsno , & ! volume per unit area of snow (m) - zlvl , & ! atm level height (m) + zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) uatm , & ! wind velocity components (m/s) vatm , & wind , & ! wind speed (m/s) @@ -2177,7 +2179,8 @@ icepack_step_therm1 real (kind=dbl_kind), optional, intent(in) :: & HDO_ocn , & ! ocean concentration of HDO (kg/kg) H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) - H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) + H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) + zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice