diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90 index a84e7d967..06012eb63 100644 --- a/columnphysics/icepack_algae.F90 +++ b/columnphysics/icepack_algae.F90 @@ -67,6 +67,9 @@ module icepack_algae private public :: zbio, sklbio + real (kind=dbl_kind), parameter :: & + exp_argmax = c10 ! maximum argument of exponential + !======================================================================= contains @@ -903,8 +906,8 @@ subroutine z_biogeochemistry (n_cat, dt, & Sink_top, & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s) ocean_b, & ! ocean_bio sum_react, & - rtau_ret, & ! retention frequency (s^-1) - rtau_rel , & ! release frequency (s^-1) + exp_ret, & ! exp dt/retention frequency + exp_rel, & ! exp dt/release frequency atm_add_cons , & ! zbgc_snow+zbgc_atm (mmol/m^3*m) dust_Fe , & ! contribution of dust surface flux to dFe (umol/m*3*m) source , & ! mmol/m^2 surface input from snow/atmosphere @@ -933,7 +936,8 @@ subroutine z_biogeochemistry (n_cat, dt, & ! when hin > hbri: just used in sw calculation real (kind=dbl_kind):: & - bio_tmp ! temporary variable + bio_tmp, & ! temporary variable + exp_min ! temporary exp var real (kind=dbl_kind):: & Sat_conc , & ! adsorbing saturation concentration (mmols/m^3) @@ -1058,21 +1062,27 @@ subroutine z_biogeochemistry (n_cat, dt, & !----------------------------------------------------------------- ! time constants for mobile/stationary phase changes !----------------------------------------------------------------- - + + exp_rel(m) = c0 + exp_ret(m) = c0 + if (tau_ret(m) > c0) then + exp_min = min(dt/tau_ret(m),exp_argmax) + exp_ret(m) = exp(-exp_min) + endif + if (tau_rel(m) > c0) then + exp_min = min(dt/tau_rel(m),exp_argmax) + exp_rel(m) = exp(-exp_min) + endif if (m .ne. nlt_bgc_N(1)) then if (hin_old > hin) then !melting - rtau_rel(m) = c1/tau_rel(m) - rtau_ret(m) = c0 + exp_ret(m) = c1 else !not melting - rtau_ret(m) = c1/tau_ret(m) - rtau_rel(m) = c0 + exp_rel(m) = c1 endif elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then - rtau_rel(m) = c1/tau_rel(m) - rtau_ret(m) = c0 + exp_ret(m) = c1 elseif (tr_bgc_N) then - rtau_ret(m) = c1/tau_ret(m) - rtau_rel(m) = c0 + exp_rel(m) = c1 endif ocean_b(m) = ocean_bio(m) @@ -1152,8 +1162,8 @@ subroutine z_biogeochemistry (n_cat, dt, & do k = 1,nblyr+1 initcons_mobile(k) = in_init_cons(k,mm)*trcrn(nt_zbgc_frac+mm-1) initcons_stationary(k) = mobile(mm)*(in_init_cons(k,mm)-initcons_mobile(k)) - dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp(-dt*rtau_ret( mm))-c1) + & - initcons_stationary(k)*(c1-exp(-dt*rtau_rel(mm)))) + dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp_ret(mm)-c1) + & + initcons_stationary(k)*(c1-exp_rel(mm))) initcons_mobile(k) = max(c0,initcons_mobile(k) + dmobile(k)) initcons_stationary(k) = max(c0,initcons_stationary(k) - dmobile(k)) if (initcons_stationary(k)/hbri_old > Sat_conc) then @@ -1270,7 +1280,7 @@ subroutine z_biogeochemistry (n_cat, dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'in_init_cons(:,mm):',in_init_cons(:,mm) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'rtau_ret( mm),rtau_rel( mm)',rtau_ret( mm),rtau_rel( mm) + write(warnstr,*) subname, 'exp_ret( mm),exp_rel( mm)',exp_ret( mm),exp_rel( mm) call icepack_warnings_add(warnstr) write(warnstr,*) subname,'darcyV,dhtop,dhbot' call icepack_warnings_add(warnstr) @@ -1370,7 +1380,7 @@ subroutine z_biogeochemistry (n_cat, dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'rtau_ret( m),rtau_ret( m)',rtau_ret( m),rtau_ret( m) + write(warnstr,*) subname, 'exp_ret( m),exp_ret( m)',exp_ret( m),exp_ret( m) call icepack_warnings_add(warnstr) call icepack_warnings_add(subname//'negative bgc') call icepack_warnings_setabort(.true.,__FILE__,__LINE__) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index f26f905e3..77e73719f 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -161,7 +161,7 @@ subroutine atmo_boundary_layer (sfctype, & psimhu, & ! unstable part of psimh psixhu ! unstable part of psimx - character(len=*),parameter :: subname='(atm_boundary_layer)' + character(len=*),parameter :: subname='(atmo_boundary_layer)' !------------------------------------------------------------ ! Define functions @@ -389,8 +389,7 @@ subroutine atmo_boundary_const (sfctype, calc_strair, & Tsf, potT, & Qa, & delt, delq, & - lhcoef, shcoef, & - Cdn_atm) + lhcoef, shcoef ) character (len=3), intent(in) :: & sfctype ! ice or ocean @@ -407,9 +406,6 @@ subroutine atmo_boundary_const (sfctype, calc_strair, & wind , & ! wind speed (m/s) rhoa ! air density (kg/m^3) - real (kind=dbl_kind), intent(in) :: & - Cdn_atm ! neutral drag coefficient - real (kind=dbl_kind), intent(inout):: & strx , & ! x surface stress (N) stry ! y surface stress (N) @@ -429,7 +425,7 @@ subroutine atmo_boundary_const (sfctype, calc_strair, & tau, & ! stress at zlvl Lheat ! Lvap or Lsub, depending on surface type - character(len=*),parameter :: subname='(atm_boundary_const)' + character(len=*),parameter :: subname='(atmo_boundary_const)' !------------------------------------------------------------ ! Initialize @@ -901,6 +897,8 @@ subroutine icepack_atm_boundary(sfctype, & worku = uvel endif + Cdn_atm_ratio_n = c1 + if (trim(atmbndy) == 'constant') then call atmo_boundary_const (sfctype, calc_strair, & uatm, vatm, & @@ -909,8 +907,7 @@ subroutine icepack_atm_boundary(sfctype, & Tsf, potT, & Qa, & delt, delq, & - lhcoef, shcoef, & - Cdn_atm) + lhcoef, shcoef ) if (icepack_warnings_aborted(subname)) return else ! default call atmo_boundary_layer (sfctype, & diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90 index ea6847b1f..657b9965a 100644 --- a/columnphysics/icepack_brine.F90 +++ b/columnphysics/icepack_brine.F90 @@ -8,7 +8,7 @@ module icepack_brine use icepack_kinds - use icepack_parameters, only: p01, p001, p5, c0, c1, c2, c1p5, puny + use icepack_parameters, only: p01, p001, p5, c0, c1, c2, c1p5, c10, puny use icepack_parameters, only: gravit, rhoi, rhow, rhos, depressT use icepack_parameters, only: salt_loss, min_salin, rhosi use icepack_parameters, only: dts_b, l_sk @@ -45,6 +45,9 @@ module icepack_brine b1 = 1000.0_dbl_kind, & ! (kg/m^3) b2 = 0.8_dbl_kind ! (kg/m^3/ppt) + real (kind=dbl_kind), parameter :: & + exp_argmax = 30.0_dbl_kind ! maximum argument of exponential for underflow + !======================================================================= contains @@ -511,6 +514,7 @@ subroutine update_hbrine (meltb, meltt, & darcy_coeff, & ! magnitude of the Darcy velocity/hbrocn (1/s) hbrocn_new , & ! hbrocn after flushing dhflood , & ! surface flooding by ocean + exp_arg , & ! temporary exp value dhrunoff ! direct runoff to ocean real (kind=dbl_kind), parameter :: & @@ -540,7 +544,13 @@ subroutine update_hbrine (meltb, meltt, & bphi_min = bphin dhrunoff = -dhS_top*aice0 hbrocn = max(c0,hbrocn - dhrunoff) - hbrocn_new = hbrocn*exp(-darcy_coeff/bphi_min*dt) + exp_arg = darcy_coeff/bphi_min*dt +! tcx tcraig avoids underflows but is not bit-for-bit +! if (exp_arg > exp_argmax) then +! hbrocn_new = c0 +! else + hbrocn_new = hbrocn*exp(-exp_arg) +! endif hbr = max(hbrmin, h_ocn + hbrocn_new) hbrocn_new = hbr-h_ocn darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index d36bbc43a..4087b81df 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -678,7 +678,7 @@ subroutine init_vertical_profile(nilyr, nslyr, & zSin ! internal ice layer salinities real (kind=dbl_kind), dimension (:), & - intent(out) :: & + intent(inout) :: & zqsn , & ! snow enthalpy zTsn ! snow temperature diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index 7c2194ca7..5b480de96 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -135,7 +135,7 @@ subroutine init_restart use icedrv_init_column, only: init_hbrine, init_bgc use icedrv_restart, only: restartfile use icedrv_restart_shared, only: restart - use icedrv_restart_column, only: read_restart_bgc + use icedrv_restart_bgc, only: read_restart_bgc use icedrv_state ! almost everything integer(kind=int_kind) :: & diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index ba774d166..3507444a5 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -94,7 +94,7 @@ subroutine ice_step use icedrv_flux, only: scale_factor, init_history_therm, init_history_bgc, & daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd, init_history_dyn use icedrv_restart, only: dumpfile, final_restart - use icedrv_restart_column, only: write_restart_bgc + use icedrv_restart_bgc, only: write_restart_bgc use icedrv_state, only: trcrn use icedrv_step, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_ridge, step_radiation, & diff --git a/configuration/driver/icedrv_restart_column.F90 b/configuration/driver/icedrv_restart_bgc.F90 similarity index 99% rename from configuration/driver/icedrv_restart_column.F90 rename to configuration/driver/icedrv_restart_bgc.F90 index 26ccdb552..807ead8bc 100644 --- a/configuration/driver/icedrv_restart_column.F90 +++ b/configuration/driver/icedrv_restart_bgc.F90 @@ -4,7 +4,7 @@ ! ! author: Elizabeth C. Hunke, LANL ! - module icedrv_restart_column + module icedrv_restart_bgc use icedrv_kinds use icedrv_constants @@ -20,7 +20,8 @@ module icedrv_restart_column implicit none private - public :: write_restart_bgc, read_restart_bgc + public :: write_restart_bgc, & + read_restart_bgc !======================================================================= @@ -667,6 +668,6 @@ end subroutine read_restart_bgc !======================================================================= - end module icedrv_restart_column + end module icedrv_restart_bgc !======================================================================= diff --git a/configuration/scripts/tests/report_results.csh b/configuration/scripts/tests/report_results.csh index cf18261a8..ad37f2eb4 100755 --- a/configuration/scripts/tests/report_results.csh +++ b/configuration/scripts/tests/report_results.csh @@ -148,8 +148,8 @@ if ( ${case} =~ *_${compiler}_* ) then if (${fcomp} == "") set rcomp = ${gray} if (${ftime} == "") set rtime = ${gray} - set fcomp = `grep " ${case} " results.log | grep " bfbcomp" | grep "baseline-does-not-exist" | wc -l ` - if ($fcomp > 1) set rcomp = ${gray} + set fregrx = `grep " ${case} " results.log | grep " compare" | grep "baseline-does-not-exist" | wc -l ` + if ($fregrx > 0) set rregr = ${gray} if (${rbuild} == ${red}) set tchkpass = 0 if (${rrun} == ${red}) set tchkpass = 0 diff --git a/doc/source/developer_guide/dg_driver.rst b/doc/source/developer_guide/dg_driver.rst index 945cacd1f..b8389b90f 100755 --- a/doc/source/developer_guide/dg_driver.rst +++ b/doc/source/developer_guide/dg_driver.rst @@ -29,7 +29,7 @@ The icepack driver consists of the following files | **icedrv_init.F90** general initialization routines | **icedrv_init_column.F90** initialization routines specific to the column physics | **icedrv_restart.F90** driver for reading/writing restart files -| **icedrv_restart_column.F90** restart routines specific to the column physics +| **icedrv_restart_bgc.F90** restart routines specific to the column physics | **icedrv_restart_shared.F90** code shared by all restart options | **icedrv_state.F90** essential arrays to describe the state of the ice | **icedrv_step.F90** routines for time stepping the major code components